Effect of Climate Change on Water Temperature and Stratification of a Small, Temperate, Karstic Lake (Lake Kozjak, Croatia)

As closed systems, lakes are extremely vulnerable to climate change. Understanding the response to climate change is crucial for effective management and conservation of the lakes and their associated ecosystems. This study focuses on Lake Kozjak, Croatia, a small lake belonging to the Plitvice Lakes system. This system represents a unique hydrogeological karstic phenomenon, closely dependent on a delicate biochemical balance necessary for tufa formation. We apply a simple one-dimensional model, SIMO v.1.0, to predict future water temperature in Lake Kozjak under three scenarios (RCP2.6, RCP4.5 and RCP8.5) from 2006 to 2100. The model was calibrated using measured water temperature profiles and meteorological data from a nearby station. In addition to analyzing the average temperatures of the epilimnion, hypolimnion and the whole lake, we also studied the surface and bottom layer temperatures and their relation to specific forcing parameters. The Schmidt stability index was used as a quantitative indicator to assess lake stability. The simulation results indicate average lake water temperature increase of 0.51, 1.41 and 4.51 °C (100 y)−1 for RCP2.6, RCP4.5 and RCP8.5, respectively. This increase in the water temperature is not accompanied by a substantial strengthening of stratification under RCP2.6 and RCP4.5 scenarios due to the temperature raise being present both in the epilimnion and hypolimnion. However, significant lengthening of the stratification period is observed even for the most stringent scenario, 16, 28 and 47 d (100 y)−1 for RCP2.6, RCP4.5 and RCP8.5, respectively. The predicted water temperature increase and prolonged stratification period may carry serious ecological and environmental implications. • Mean lake water temperature is projected to increase by 0.51 to 4.51 °C (100 y)−1. • Baseline scenario surface temperature increase of 5.2 °C (100 y)−1 is predicted. • Stratification period is predicted to lengthen by 16 (RCP2.6) to 47 days (RCP8.5). • Substantial stratification strengthening is expected only under RCP8.5.


Introduction
Climate change is introducing significant effects and serious threats to aquatic systems worldwide (Mooij et al. 2005;Verburg and Hecky 2009;Sahoo et al. 2013;Råman Vinnå et al. 2021).Lake responses to climate change are complex and not associated with changes in single climatic variables but a combination of factors that govern the heat budget, such as air temperature, humidity, shortwave and longwave radiation, wind speed, and precipitation.In addition, during the last couple of decades, anthropogenic pressure on aquatic systems has also globally increased (Amiri et al. 2014;Baier et al. 2014;Janse et al. 2015;Gude 2017).To better understand and predict possible lake responses and eventual consequences, increasing effort is being invested into researching not only the surface water response but also the deep water response and stratification, including studies dealing with specific lakes (Matsumoto et al. 2019;Ayala et al. 2020;Piccioni et al. 2021;Noori et al. 2022b;Sterckx et al. 2023), but also studies with regional or global extent (Woolway and Merchant 2019;Shatwell et al. 2019;Maberly et al. 2020;Pilla et al. 2020;Råman Vinnå et al. 2021;Woolway et al. 2021Woolway et al. , 2022;;Golub et al. 2022;Sharaf et al. 2023;Woolway 2023).While surface water temperatures have exhibited a globally increasing trend (O'Reilly et al. 2015), deep water temperatures do not follow this consistent trend and show high variability in both sign and value across lakes (Pilla et al. 2020).Read et al. (2014) also stress the contrast between the strongly coherent surface temperatures opposed to weakly coherent hypolimnion temperatures and stratification patterns and they recognize the need to model individual lakes to capture the spatial and temporal heterogeneity in the lake responses.
Changes in lake water temperature have widespread consequences, as they directly influence biochemical processes, such as the solubility of gases and minerals and the rate of chemical reaction (Benson and Krause 1980;Krumgalz 2018;Jane et al. 2021) or methane production in bottom waters (Jansen et al. 2022).Furthermore, the strength of lake stratification and the length of the stratification period are key elements in determining the transportation and distribution of gases and nutrients within a lake (Vachon et al. 2019;Ladwig et al. 2021).Consequently, water temperature changes may affect the biodiversity in a lake (leading to eutrophication, destabilization of macrophyte-dominated clear water lakes, cyanobacteria dominance in phytoplankton communities, shifts in habitats, promotion of invasive species, etc.) and its environment (reducing the number of birds, affecting the mosquito populations and mosquito-borne diseases, etc.) (Mooij et al. 2005;Rasconi et al. 2017;Scholz et al. 2017;Kraemer et al. 2021;Estrela-Segrelles et al. 2023).
This study aims to contribute to global research on lakes by analyzing the response of a small, karstic lake (Lake Kozjak) to climate change.This lake is located in the mountainous region of Croatia and is one of the karstic Plitvice Lakes, which is a system of lakes that are interconnected by cascades and waterfalls.The lake chain consists of 16 named and several unnamed cascading lakes, with a total water surface of 1.949 km 2 .The first, southernmost lake in the chain (Prošće Lake) is at 636 m above sea level (a.s.l.) while the last lake (Novakovića Brod Lake) is at 503 m a.s.l.The lake system stretches in an approximate south-north direction.The territory of Plitvice Lakes National Park is part of the Dinaric karst area.It encompasses highly distinctive geological, geomorphological and hydrological features that create the extraordinary beauty of the landscapes.The carbonate rocks dominating the terrain were vital in the formation of the specific and extremely rich morphology in the area.In the park area, there are approximately 8000 recorded sinkholes, several karstic plateaus and 114 speleological sites with total explored length of 1664 m and depth of 2251 m (Kovačević et al. 2019).The carbonate rocks also have a significant influence on the mineral content of the water (Petalas et al. 2018;Pratama et al. 2021).The calcium carbonate dissolved in the water is one of the primary conditions allowing the tufa formation process responsible for the formation of lakes, cascades and waterfalls.Thanks to their natural beauty, exceptional biodiversity and hydrogeological uniqueness, the Plitvice Lakes were the first National Park in Croatia, established in 1949(PLNP 2023a) and have been inscribed on the UNESCO World Heritage List since 1979 (UNESCO 2023).Their uniqueness is primarily due to the biodynamic processes of tufa formation, which are responsible for the ongoing creation of natural dams and barriers that shape the lake system.These processes depend on a fragile balance of physical, chemical and biological conditions.They occur in a narrow range of water temperatures and chemical conditions in the presence of certain living organisms (i.e., microphytes, macrophytes and bacteria).Other important factors include the quantity and velocity of water flow.
In this study, the response of the largest of the Plitvice Lakes, Lake Kozjak, to future climate changes is evaluated for the first time.For the prediction of the vertical temperature profiles, we use a simple 1-D energy budget model (SIMO) we had previously developed and validated against lake temperature data measured continuously during an observational campaign in two lakes belonging to the Plitvice Lakes, Lake Kozjak and Lake Prošće (Šarović et al. 2022).We project the lake water temperature and stratification under three different Representative Concentration Pathway (RCP) scenarios.The RCP scenarios are greenhouse concentration scenarios used for climate research in the IPCC fifth Assessment Report (IPCC 2014).There are four scenarios named according to the radiative forcing in 2100 (2.6, 4.5, 6, and 8.5 W m −2 ).RCP2.6 is the most stringent pathway that requires negative CO 2 emissions and is likely to maintain the global temperature rise below 2 °C by 2100.RCP4.5 is an intermediate scenario and is the most likely.Assuming no emission reduction policies and the eventual exhaustion of nonreusable fuels, this scenario requires carbon dioxide (CO 2 ) emissions to start declining by 2045 and is expected to result in a global temperature rise between 2 and 3 °C by 2100.RCP6 assumes an emission peak in 2080.It is based on high emissions and stabilization of radiative forcing after 2100 due to different technologies and policies for emission reduction.It is expected to result in a global temperature rise of between 3 and 4 °C.RCP8.5 is the so-called "baseline scenario", which assumes no efforts to constrain emissions, which means that emissions continue to rise throughout the twenty-first century.It is expected to result in a global temperature rise of between 4 and 5 °C by 2100.RCP8.5 is not considered likely, but it is analyzed as a worst-case scenario.In this study, RCP2.6, 4.5 and 8.5 were considered.
The purpose of this study is to illuminate the potentially specific response of Lake Kozjak to climate change, taking into account its morphology (small, but relatively deep) and geographical location.The findings do not only provide valuable input information for numerous multidisciplinary studies and aid the local community in managing and protecting the lake, but they also offer valuable insights into the expected behavior of temperate lakes with similar characteristics, particularly those situated in mountainous regions.

Study Area
Lake Kozjak (Fig. 1) is the 12 th lake in the Plitvice Lakes chain, situated at 535 m a.s.l.It is the largest lake in the system, both in area (0.82 km 2 ) and volume (0.01271 km 3 ), as well as the deepest lake, with maximum and average depths of 46 m and 17.3 m, respectively.The maximum fetch of the lake is approximately 2300 m, and its width is between 100 and 600 m.A submerged barrier, whose top is 5-6 m below the surface, divides the lake into two subbasins: a deeper northern subbasin with a fetch of ~ 730 m and a shallower southern subbasin with a fetch of ~ 1570 m (Fig. 1c).The complex morphology of the lake was formed approximately 400 years ago when two smaller lakes merged due to submergence of the dam separating them.There are two theories explaining this event.The older theory (Petrik 1958;Srdoč et al. 1985) suggests that the downstream barrier (the current dam at Lake Kozjak) grew faster than the upstream barrier between the two lakes, causing the water level of the second lake to rise until it reached and surpassed the upstream barrier.This theory was replaced by a new theory in 2019, when divers discovered a point where part of the submerged barrier had probably broken off, creating a 4 m wide channel, making the downstream barrier the only active barrier and creating Lake Kozjak (PLNP 2023b).
Main source of water for Lake Kozjak is the upstream lake with an average annual flow of 2.14 m 3 s −1 (min.flow 0.6 m 3 s −1 and max.flow 11.6 m 3 s −1 ), but it also receives a substantial quantity of water from Rječica River and streams Jasenovac and Matijaševac, with total annual average flow of 0.67 m 3 s −1 (min.flow 0.2 m 3 s −1 and max.flow 3.4 m 3 s −1 ) (Biondić et al. 2010).
During the summer, the water flow through the lake system gets very low, especially in the southern part of the system, from the first lake, Lake Prošće, to Lake Kozjak (Petrik 1958).This means that the thermal processes mainly occur locally in each lake separately.The existing flow is in a form of shallow waterfalls and cascades.Therefore, during summer stratification, only water from the heated up epilimnion layer flows from the upstream to the downstream lake.Petrik (1961) also noticed that in the warm season there is a surface water temperature rise downstream of the water sources.This downstream surface water temperature rise is most pronounced in Lake Prošće and less so between Lake Prošće and Lake Kozjak exit.During the cold season, downstream temperature decrease was observed (Petrik 1961; for details refer to Klaić et al. 2018).
A number of studies using both measured historical data and predicted future scenarios report worldwide changes in lake mixing regime due to climate change (Woolway and Merchant 2019;Shatwell et al. 2019;Råman Vinnå et al. 2021).Råman Vinnå et al. (2021) underline the low to mid altitude (< 1500 m a.s.l.) small (< 0.5 km 3 ) lakes, such as Lake Kozjak, as especially sensitive to regime changes.Previous studies, generally based on older sporadic observations, classify Lake Kozjak as dimictic (Babinka 2007;Špoljar et al. 2007;Klaić et al. 2018).However, in recent years, no freezing over of the lake was observed (Plitvice Lakes National Park, personal communication).Late freezing, earlier ice loss or complete absence of ice is recognized as a widespread consequence of global warming (Sharma et al. 2019(Sharma et al. , 2021;;Imrit and Sharma 2021) with cascading effects on water temperature and stratification phenology (Woolway et al. 2021), as well as on ecological processes (Hebért et al. 2021).Moreover, continuous observation of Lake Kozjak water temperature profile from 2018 to 2020 also indicated that the lake had started to become monomictic (Šarović et al. 2022), most likely as a result of climate change.

Numerical Model
For the prediction of the vertical temperature profile in Lake Kozjak, the previously developed model, SIMO v.1.0,was used (Šarović et al. 2022).This is a simple one-dimensional model for monomictic lakes, based on the one-dimensional energy balance equation, assuming a constant horizontal cross-sectional area of the lake: where c p is the water specific heat capacity (J kg −1 K −1 ), ρ is the water density (kg m −3 ) calculated by the Chen and Millero (1986) formula assuming zero salinity, T is the water temperature (°C), t is time (s), z is depth (m), k m and k t are the molecular and turbulent thermal conductivity (W m −1 K −1 ), respectively, Φ is the heat flux (W m −2 ) and M conv is the convective mixing term (W m −3 ).The main forcing term in Eq. ( 1) is the heat source term, accounting for the meteorological forcing.It represents the sum of the net shortwave radiation (S n ), net longwave radiation (L n ), sensible heat flux (H s ), latent heat flux (H l ), and heat flux induced by precipitation (H p ).At the surface, it is calculated as: All the terms in Eq. ( 2) are completely absorbed at the lake surface layer, except for shortwave radiation, which partially passes through the water.Therefore, the heat flux at a certain depth consists of the attenuated shortwave radiation, which is calculated using the arctangent model (Henderson-Sellers 1986): where S n (z) is the net shortwave radiation at water depth z (W m −2 ), α is the water surface albedo, S is the total shortwave radiation at the lake surface (W m −2 ), and K 1 , K 2 and K 3 are empirical constants.For further details on the model and on calculating the heat flux terms, turbulent thermal conductivity, and light attenuation, please refer to Šarović et al. (2022).
In Šarović et al. (2022), the model was run without lake-specific calibration to show its general applicability.It was evaluated against measured data for two of the Plitvice Lakes (Lakes Kozjak and Prošće) and showed satisfactory performance.In this study, the model parameters were calibrated specifically for Lake Kozjak using observation data of the water temperature profile measured at the deepest point of the lake (44.8902°N, 15.6038°E) in ( 1) the period between 7 July 2018 and 2 November 2020.The water temperature was measured at 16 different depths, namely, 0.2, 0.5, 1, 1. 5, 3, 5, 7, 9, 11, 13, 15, 17, 20, 23, 27 and 43 m.The data were measured and logged with HOBO TidBiT MX Temp 400, as previously described in Klaić et al. (2020).The accuracy of the sensors is ± 0.20 °C for temperatures between 0 and 70 °C and ± 0.25 °C for temperatures between − 20 and 0 °C.The sampling frequency was 1 Hz and 2 min means were stored.For the purpose of this study hourly mean lake temperatures were used.
Hourly meteorological input data used to force the calibration simulation were obtained from the Plitvička Jezera automatic meteorological station (44.8811°N, 15.6197° E), located approximately 200 m from the Lake Kozjak shoreline and 1.6 km southeast of the water temperature measuring site.The station is part of the Croatian Meteorological and Hydrological Service (CMHS) network.The meteorological and water temperature data availability periods do not completely overlap due to gaps in the meteorological data (4 November 2018 -1 January 2019, 31 December 2019 -2 July 2020).Taking into account the data availability, the calibration was performed based of a year-long simulation with hourly time step, initiated on 1 January 2019 with a measured water temperature profile.
Since the solar heat flux was identified as the main factor affecting the lake water temperature, the main focus was on the adjustment of the net shortwave radiation flux and light extinction coefficient parameters.After the calibration, the light extinction coefficient was increased from 0.1 (generic value appropriate for clear oligotrophic lakes used in Šarović et al. 2022) to 0.15.Average Secchi depth of the lake, measured monthly during 2009 and 2010, was z SD = 10.5 m (Žutinić et al. 2014) which indicates a light extinction coefficient value of 1.7 /z SD = 0.16 (Poole and Atkins 1929).The calibrated light extinction value roughly agrees with the measured data from 2009 -2010.
As an additional improvement, the lake cross-sectional area shrinkage with depth was taken into account in the convective mixing algorithm.Since precise bathymetry data was not available, it was assumed that the cross-sectional area decreases linearly with depth.This successfully removed the pronounced temperature overprediction during winter convective overturn (Figs. 2 and 3).On the other hand, it introduced an underestimation of the temperature in the topmost layers, where the cross-sectiona does not follow the linear assumption, during the short and shallow convective mixing occasionally happening in January and February.This underestimation is present only in the topmost layers during a relatively short periods of time.Moreover, these kind of conditions are not occurring often in the analyzed future scenarios.Šarović et al. (2022) reported bias and root mean square error (RMSE) of uncalibrated SIMO, calculated from a year-long simulation, of 1.4 and 1.9 °C, respectively.After the performed calibration these values were reduced to 0.1 and 0.7 °C, respectively.
Additionally, in Šarović et al. (2022), the Lake Kozjak model was discretized using 16 layers corresponding to the measuring depths of the available data, with a vertical resolution that decreased with depth.In this study, the resolution was kept at a minimum of 2 m, even in the deeper layers, to more precisely capture the stratification process and deepening of the thermocline.In total, 25 layers were used.The depths of the integration points were 0.2, 0.5, 1, 1.5 and 3 m, and after the 3 m depth they were spaced at 2 m intervals.A temporal resolution of 1 h was used.
Note that the lake has been modeled as isolated from the lakes in the chain.The satisfactory model performance presented in Šarović et al. (2022) proved that this simplification is reasonable.

Forcing Data
Projected meteorological data for the period between 2006 and 2100 were obtained from the Coordinated Regional Climate Downscaling Experiment (CORDEX).CORDEX is an international coordinated framework initiated by the World Climate Research Programme (WCRP) as an effort to evaluate and improve regional climate downscaling (RCD) techniques and produce a new generation of RCD-based fine-scale climate projections (COR-DEX 2022).
In the present study, air temperature, wind speed, air humidity, atmospheric pressure, precipitation (Fig. 4a-e) and shortwave radiation data from simulations of the European domain with spatial resolutions of 0.11° and temporal resolutions of 1 h were downloaded from the Earth System Grid Federation (ESGF) node hosted by the German Climate Computing Centre (ESFG 2022).The data were from the r1i1p1 ensemble and v1 downscaling realization, obtained by downscaling the output of the NCC NorESM1-M general circulation model (GCM) for RCP2.6 and RCP4.5, and the IPSL CM5A-MR GCM for RCP8.5, using the REMO2015 regional climate model (RCM).The values for the exact investigated location were calculated via bilinear interpolation using the data from the four closest available points.An altitude difference correction assuming adiabatic lapse was applied to the air temperature values before interpolation.The shortwave radiation data were used for calculating the daily radiation profile only, while the actual total daily shortwave irradiation was calculated by SIMO v.1.0as described in previous research (Šarović et al. 2022).The calculated shortwave radiation used as an input in the water temperature prediction model is shown in Fig. 4f.The climate projections involve inherent uncertainties originating from multiple sources, usually partitioned as internal, scenario and model uncertainty.Internal uncertainty is present due to internal unforced variability, it is inherent to the climate system and cannot be reduced.Scenario uncertainties are related to factors affecting future emissions (population, technology, socio-political situation, etc.).Model uncertainties are consequence of the simplifications, parametrizations and approximations used in the model.Model uncertainties are introduced in two steps, when running the GCM and later during downscaling with RCM.Since the choice of GCM generally has larger influence on the climate projections than the choice of RCM (Christensen and Kjellström 2020), it would have been more appropriate for this study to use results obtained by same GCM for all three scenarios.However, such data set, containing all the necessary variables with the desired spatial and temporal resolution, was not available.With respect to the temperature projections, in the area of Lake Kozjak, the NorESM1-M downscaled with REMO2015 is positively biased (~ 2 °C), but has negatively biased temperature response to climate change (~ -1 °C), Fig. 3 Hourly values of the predicted water temperature for Lake Kozjak in 2019 at depths corresponding to the temperature sensor depths before (blue) and after (red) model calibration compared to the observed values (black) (Christensen and Kjellström 2020;Vautard et al. 2021).The opposite is true for CM5A-MR downscaled with REMO2015 -it is negatively biased (~ -2 °C), but has positively biased temperature response to climate change (~ 1 °C).Both models tend to underpredict precipitation in the winter (bias ~ -1.5 mm d −1 ) and overpredict it in the summer season (bias ~ -1 mm d −1 ).Both models generally underpredict wind speed (bias ~ -0.4 m s −1 ).On the other hand, the climate projections probably overestimate the wind speed at the given location as they do not include the fine scale orography surrounding and sheltering Lake Kozjak.Namely, wind speeds measured in the period between July 2018 and September 2020 were mainly below 2 m s −1 (Šarović et al. 2022).

Lake Water Temperature and Stratification Trend Estimations
Annually averaged water temperature trends were analyzed at different depths.From the temperature profiles, the average temperatures of the epilimnion, hypolimnion and whole lake were calculated and analyzed.The thermocline depth was determined as the depth of the highest vertical water temperature gradient, considering depths over 1 m to exclude high gradients that may be encountered in the top layers.In calculating the yearly averaged thermocline depth, only values obtained during periods of lake stratification were taken into account.The first and last layer in the discretized mode, ranging from 0 to 0.35 m and from 42 to 46 m depth, respectively, represent the surface and bottom layer.Their water temperature evolution was more closely examined by using the monthly averaged water temperature and calculating the temperature change rate for each month separately.In this way, the climate change effect on the yearly variation was examined.
The Schmidt stability index (SSI) was used as a quantitative indicator of lake stability.This index represents the amount of energy required to mix the water column to obtain a uniform water temperature profile without heat exchange.Assuming a constant lake crosssectional area, the SSI was calculated using the formula presented in Gaudard et al. (2019): where g is the acceleration of gravity (m s −2 ), z denotes the depth (m), z is the mean lake depth (m), and ρ z is the water density at depth z (kg m −3 ).The SSI was also used to determine the beginning, end and duration of the stratification period.The timing of the beginning and end of stratification was expressed in day of the year.Gaudard et al. (2019) suggested SSI/z lake ≥ 10 J m −3 as a stratification criterion.For Lake Kozjak, this criterion translates to an SSI ≥ 460 J m −2 .
Note that the change rates of all relevant variables are calculated over the whole simulation period, from 2006 to 2100, and they are expressed as change per 100 years (e.g., °C (100 y) −1 ).

Results
Figure 5 shows the yearly averaged values of the water temperature and the vertical gradient at different depths for the three different scenarios, namely, RCP2.6,RCP4.5, and RCP8.5.While there was an obvious increase in the water temperature under scenario RCP4.5 (Fig. 5b1) and even more so under RCP8.5 (Fig. 5c1), it was not accompanied by a significant increase in the vertical temperature gradient (Fig. 5a2-c2) as the temperature increase occurred at all depths.Furthermore, there is no significant trend of the depth of the maximum vertical gradient indicating no significant trend of the average thermocline depth.As it will be later demonstrated in more detail, significant increase in the maximum thermocline depth is observed under scenario RCP8.5, along with extension of the stratification period in all three analyzed scenarios.However, these findings are not reflected in Fig. 5 due to the annual averaging.
To further examine the stratification evolution, linear regression of the yearly averaged water temperature was calculated for each depth and scenario.Figure 6a1-c1 show the yearly averaged water temperature in each layer, along with the thermocline temperature.The thermocline temperature is defined as the temperature at the depth of the highest temperature gradient.Since the temperature gradient is calculated at the border between two layers, the thermocline temperature is calculated as the temperature average of those two layers.Since Fig. 6a1-c1 show yearly averaged temperature values, the thermocline temperature was also calculated as a yearly average using the whole year data and not just data from the periods when the lake was stratified which is generally not the usual approach; so, this should be taken into account if the results are compared to other studies.The thermocline temperature exhibits a significant increase trend of 0.8, 1.82 and 4.68 °C (100 y) −1 under scenarios RCP2.6,RCP4.5 and RCP8.5, respectively.
The temperature change rate coefficients were also calculated for each layer and they are shown as a function of the water depth in Fig. 6a2 (RCP8.5), the water temperature increased at a much higher rate than in the other two scenarios.This was expected considering the significantly higher air temperature and shortwave radiation increase rate and relative humidity decrease rate in that scenario compared to the other two scenarios (Fig. 4a, c and f).Surface layer water temperature increase rates of 0.62, 1.75 and 5.19 °C (100 y) −1 and bottom layer water temperature increase rates of 0.22, 1.05 and 4.06 °C (100 y) −1 were calculated for scenarios RCP2.6,RCP4.5 and RCP8.5, respectively.However, trend significance check showed that the calculated trends for depths greater than 20 m under scenario RCP2.6 (Fig. 6a2) were not significant (p > 0.05).Generally, the water temperature change rate in the hypolimnion decreased with depth in a roughly linear manner (Fig. 6a2-c2).This was not the case in the epilimnion and the thermocline region due to the effect of different meteorological forcing mechanisms and the nonlinearity of the processes governing the lake thermal cycle (thermal stratification and convective overturn).In the shallow layers (depths < 10 m), the effect of the change in shortwave radiation could be observed as a significantly higher change rate under scenario RCP8.5 and a slower change rate under scenario RCP2.6 compared to those in the slightly deeper layers (depths 10 to 15 m).As no significant thermocline deepening was expected under scenario RCP2.6, there was a steep decrease in the change rate through the epilimnion and thermocline region (up to ~ 23 m); and in the hypolimnion, the change rate only slightly decreased with depth.A steep decrease in the change rate in the epilimnion was also observed under scenarios RCP4.5 and RCP8.5; however, in these two scenarios, a sudden increase was exhibited at a depth of ~ 21 m due to the deepening of the epilimnion.
The slower temperature increase at greater depths was in agreement with previous studies (Pilla et al. 2020;Noori et al. 2022a); however, we are not aware of any studies reporting the rate as a function of water depth with a resolution that is fine enough to discern this kind of peculiar behavior in the thermocline region.
Results from a study investigating the climate change impact on 29 Swiss lakes predicted an average increase in surface water temperature of 3.3, 1.7 and 0.9 °C and bottom temperature of 1.6, 0.93 and 0.48 °C by the end of the century for RCP8.5, RCP4.5 and RCP2.6, respectively (Råman Vinnå et al. 2021).A worldwide study including 635 lakes (Woolway and Merchant 2019) predicts an average lake surface warming of 2.5 and 1.1 °C by the end of the century, for RCP6.0 and RCP2.6, respectively.However, the most extreme predicted lake surface warming by the end of the century under the scenario RCP6.0 reaches around 5.5 °C.The top layer water temperature results for Lake Kozjak are in agreement with the results from the existing studies.On the other hand, bottom layer temperatures in Lake Kozjak exceed the predictions presented in Råman Vinnå et al. (2021).The authors of that study explain the mild heating of the deep waters by the fact that the winter cooling phase remains strong enough to remove the heat accumulated in summer, and cool the entire lake down to the temperature of maximum density of ~ 4 °C, thus resetting the lake bottom temperature every winter.This is not the case in Lake Kozjak where the bottom layer temperature rarely has dropped to ~ 4 °C after 2030.
Figure 7 shows the yearly values of the depth-averaged temperature for the epilimnion, hypolimnion and whole lake.The lake average temperature was closer to the hypolimnion average temperature, as the hypolimnion was roughly twice as thick as the epilimnion.The epilimnion temperature was, as expected, significantly higher.This was also true for the linear regression coefficients.Namely, the epilimnion regression coefficients were similar to those of the air temperature (Fig. 4a).
Generally, water temperature shows lower yearly variability with depth.The bottom layer (from 42 to 46 m depth) water temperature in Lake Kozjak was close to constant throughout the year (the standard deviation calculated from measured data in 2019 was 0.4 °C).Consequently, the bottom layer temperature increase rate calculated for each month was also fairly constant (Fig. 8b).On the other hand, the surface layer (from 0 to 0.35 m depth) water temperature showed significant variability throughout the year (the standard deviation calculated from measured data in 2019 was 6.8 °C), which resulted in variability in the monthly temperature increase rate (Fig. 8a).The highest surface layer water temperature increase rate was calculated for the RCP8.5 scenario during summer stratification, reaching a maximum of 6.63 °C (100 y) −1 in August.The highest peaks for RCP2.6 and RCP4.5 occurred in April and February, respectively.The surface layer water temperature increase rate peaks corresponded to the shortwave radiation flux change rate peaks (Fig. 9c), which was the dominant component of the total heat flux during the summer stratification period (Fig. 10a).A significant part of the shortwave radiation flux (40%) is absorbed in the topmost 0.6 m of depth.Exponential decay starts below this surface absorption layer, reducing the shortwave radiation flux at 20 m to ~ 3.5% of the surface shortwave radiation flux.Obviously, the direct solar radiation effect at depths around or greater than 20 m is negligible.All the other components of the total surface heat flux except the shortwave radiation are absorbed in the surface layer.This means that the main heat transfer to the deep layers is through diffusion which is relatively slow and does not show high variation throughout the year.However, there is slightly higher temperature increase rate in the bottom layer in the winter months compared to the rest of the year and the reason for this is the winter convective overturn.
The yearly total heat flux had a positive trend of 0.02, 1.38 and 3.74 W m −2 (100 y) −1 for RCP2.6,RCP4.5 and RCP8.5, respectively (Fig. 10b), leading to an obvious increase in the mean lake water temperature.However, the total heat flux change rate that was calculated separately for different months was variable, and in some scenarios, it was even negative for some months (Fig. 9b).
Figure 9b-f show the change rate of the monthly values of the heat flux and heat flux components, and Fig. 9a shows the change rate of the monthly values of air temperature, as one of the main factors affecting the heat fluxes.The shortwave radiation change was mainly due to cloud cover and relative humidity change.The decrease in relative humidity, together with the increase in air temperature, which was especially significant in the RCP8.5 scenario, led to an increase in the outgoing latent heat flux from the lake surface (Fig. 9f).The net longwave radiation and sensible heat flux change rates were less significant (Fig. 9d and e).
As mentioned in Section 2.1, recent data suggest that Lake Kozjak has undergone a change in mixing regime, transitioning from dimictic to monomictic.In all three simulated scenarios of the present study, Lake Kozjak continued to behave as a monomictic lake -it experienced stratification during the summer period and convective overturn in autumn and winter without undergoing winter stratification and spring convective overturn.A qualitative comparison of the water temperature contour plots revealed similar thermal behavior both at the beginning and end of the analyzed period (Fig. 11a-f) and no change in the mixing regime.
The yearly maximum Schmidt stability index value showed a significant positive trend for RCP4.5 and RCP8.5 of 800 and 3398 J m −2 (100 y) −1 , respectively.The trend of the yearly averaged SSI was found significant in all three scenarios, but relatively low under scenarios RCP2.6 and RCP4.5 where its value was 176 J m −2 (100 y) −1 and 377 J m −2 (100 y) −1 , respectively.As expected, under scenario RCP8.5 the trend of the yearly average SSI is higher and equals 1270 J m −2 (100 y) −1 (Fig. 12a).On the other hand, no significant trend was detected in the yearly average depth of the thermocline (Fig. 12b).However, the maximum thermocline depth increases at a rate of 7.3 m (100 y) −1 under scenario RCP8.5.Under scenario RCP4.5, a significant trend of 3 m (100 y) −1 was detected only for the period until 2085.No significant trend was detected under scenario RCP2.6.
The increase in the SSI and maximum yearly thermocline depth (Fig. 12), as well as the increase rate of the temperature difference between the surface and the bottom layers (Fig. 8c), indicated an increase in water column stability.On the other hand, the relatively low increase in the yearly average SSI led to the conclusion that this increase in the water column stability was substantial only under scenario RCP8.5.The reason that more significant increase in stability was not observed was because of the warming up of even the deepest layers of the lake, partly due to its relatively low depth which allows for enough diffusive heat transfer towards the deeper layers and also convective mixing of the entire water column in winter.This is not the case in significantly deeper lakes where the convective overturn happens only until certain depth (e.g., Perroud et al. 2009).This was easily observed, especially for RCP8.5, by comparing Fig. 11c and f, where the increase in epilimnion temperature is clearly depicted, as is the significant warming of the hypolimnion.
The stratification criterion suggested in Section 2.4 was tested and found to be adequate for Lake Kozjak, as it satisfactorily captured the beginning of stratification in spring and the convective overturn in autumn (Fig. 11).Since the SSI represents the amount of energy required to mix the water column, the surface to bottom temperature difference corresponding to the SSI value of 460 J m −2 depends on the shape of the temperature profile.In the simulation results, the surface to bottom temperature difference values corresponding to the SSI value of 460 J m −2 range between 3.3 (P 5 ) and 6.6 °C (P 95 ), with a mean of 4.7 °C.Due to the steepness of the Schmidt stability index curve (Fig. 11g-i), the calculated values of the beginning, ending and duration of stratification were relatively insensitive to the selection of the Schmidt stability criterion value.The results are shown in Fig. 13.
Surprisingly, a significant increase rate of the stratification period duration of 16.1 d (100 y) −1 was predicted even for RCP2.6.Under this scenario, stratification is expected to start approximately 10.8 days earlier and end approximately 5.2 days later per century  (Fig. 13a).Under the RCP4.5 scenario, stratification is expected to start approximately 13.2 days earlier and end approximately 14.5 days later per century, resulting in an increase in the stratification period duration of 27.7 days (Fig. 13b).The increase in the stratification period duration is highest for the RCP8.5 scenario, with 47.1 days, due to the 28.1-day earlier onset and 19-day later ending of stratification per century (Fig. 13c).The predictions for RCP2.6 agree with the results presented in Woolway et al. (2021) where the average stratification period for Northern Hemisphere lakes is predicted to increase by 13 ± 6 days by the end of the century.On the other hand, the predictions for RCP8.5 exceed the results from the same study which predicts a stratification period extension of 33.2 ± 11.6 days by the end of the century.While the earlier onset of the stratification is in agreement with Woolway et al. (2021) where it is projected to occur 22.0 ± 7.0 days earlier, the discrepancy arises in the end of stratification which is projected to occur 11.3 ± 4.7 days later compared to the rate of 19 d (100 y) −1 predicted for Lake Kozjak.As the authors identified wind speed as one of the main factor affecting stratification break-up by enhancing turbulent heat flux towards the atmosphere and mixing the near-surface water, we assume that the reason for the bigger delay of the end of stratification in Lake Kozjak is the relatively low wind speeds in the area.

Discussion
Projected changes in lake water temperature and stratification dynamics are expected to have far reaching implications.The tufa formation, which represents a fundamental phenomenon of the Plitvice Lakes System, is positively correlated with the water temperature, but as a sensitive biodynamical process, it is also vulnerable to changes in the flow velocity and disruption of the complex balance between algae, moss and cyanobacteria (Vurnek et al. 2021;Matoničkin Kepčija and Miliša 2023).As tufa formation is faster at higher water temperatures (Srdoč et al. 1985), the results of this study, which predict a significant water temperature increase, imply that in a future, warmer climate, changes in the dynamics of tufa creation and barrier growth can be expected.This will lead to a faster evolution of and changes in the lake system morphology.
Recent studies have highlighted the role of the water temperature and stratification regime not only on the tufa forming conditions, but on the Plitvice Lakes ecosystem in  2023) deals with the water chemistry of the Plitvice Lakes and shows how water temperature and stratification affect the solubility of calcium carbonate, the concentration of dissolved oxygen, the rate of diffusion of CO 2 , the water pH and hardness.Ternjej et al. (2023) mentions that the dominance of certain phytoplankton species in Plitvice Lakes is highly sensitive to water temperature and stratification.Furthermore, when the lake is mixed, the production of chlorophyll is going on through the entire water column, while the production is confined to the epilimnion during the stratification period.A significant correlation between the total macrozooplankton abundance and water temperature was found as temperature affects growth, determines the duration of each phase in the life cycle and the occurrence of individual species throughout the year (Ternjej et al. 2023).Animals are also significantly affected by water temperature.Emergence phenology of different species of insects closely depend on the surface water temperature (Čmrlec et al. 2013;Ivković and Pont 2016;Ivković et al. 2023).Benthic macroinvertebrates enter a passive drift due to environmental changes (Miliša et al. 2023).Some fish species, such as the trout, which is one of the only three species native to the Plitvice Lakes area, require relatively cold, oxygen-rich water which (Buj et al. 2023).
Lake Prošće, the first lake in the chain, is characterized as mesotrophic, while the rest of the Plitvice Lakes are characterized as oligotrophic.However, a progress of marshland communities and aquatic vegetation related to eutrophication of oligotrophic karstic lakes has already been observed (Alegro et al. 2023).Moreover, although the Plitvice Lakes area is designated as a protected area and substantial efforts are being made to minimize human impact, the risk of cultural eutrophication persists and has not yet been fully eradicated.Intense tourism development in the Plitvice Lakes area has already been identified as a threat leading to increases in pollution and waste production (Roglić 1951;Petrik 1961;Böhm 1997), possibly leading to lake eutrophication.The unresolved sewage system as well as the unsatisfactorily resolved drainage of wastewater (affluent) from farms in the catchment area pose a main concern.Water intake with growing demand, with main structures at Lake Kozjak, also exerts significant pressure on the ecosystem.The situation is worsened by a trend of building new accommodation capacities without ensured water supply and sewage systems.A major state road also used for transport of hazardous and heavy cargo, passes through the catchment area and represents an additional pollution risk.All of these, in conjunction with climate change and water temperature increase, may lead to algal blooms and eutrophication trough promoted metabolic rate (Zhang and Chadwick 2022).The vulnerability of the system to climate and anthropogenic factors is more than obvious.
It is also important to mention the socio-economic impact of the Plitvice Lakes.Namely, the Plitvice Lakes National Park plays a crucial role in the development of the local community.It is the largest employer in the wider region, employing 30% of the permanent and 90% of the temporary employees (Kovačević et al. 2019).It also generates a substantial amount of additional jobs through various services, purchase of agricultural products and accommodation capacities.

Conclusions
The Plitvice Lakes represent a special hydrogeological karst phenomenon.Numerous studies address their physical, chemical hydrological, geological, climatological and other properties.However, to the best of our knowledge, this is the first attempt to simulate the response of one of the Plitvice Lakes to future climate scenarios.The water temperature of Lake Kozjak was simulated and analyzed under three different RCP scenarios over the period between 2006 and 2100.Increases in the mean lake water temperature of 0.51, 1.41 and 4.51 °C (100 y) −1 were calculated for the RCP2.6,RCP4.5 and RCP8.5 scenarios, respectively.The temperature increase was more significant in the epilimnion than in the hypolimnion.Analyses of the monthly averaged water temperature showed that the most significant temperature increase under the RCP8.5 scenario occurred during the summer (with a peak of 6.63 °C (100 y) −1 in August), which was not the case under the RCP2.6 and RCP4.5 scenarios.The reason for this was the significant change rate of the shortwave heat flux during the summer months under the RCP8.5 scenario.
The projected trend of the extension of the stratification period under the baseline RCP8.5 scenario was 47 d (100 y) −1 , under RCP4.5 scenario 27.7 d (100 y) −1 , and even under RCP2.6 scenario there is a significant lengthening trend of 16.1 d (100 y) −1 .No significant trend was detected in the yearly average thermocline depth, but the maximum thermocline depth value increases at a rate of 7.3 m (100 y) −1 under scenario RCP8.5 and 3 m (100 y) −1 under scenario RCP4.5.Substantial increase in the average stratification strength, evaluated based on the Schmidth stability index, was observed only in case of the baseline scenario RCP8.5.Namely, due to the relatively low depth of Lake Kozjak, there is significant warming of the deep layers through heat diffusion and convective mixing, which does not allow for more significant increase in stratification strength.
The consequences of the temperature rise and significantly longer stratification period are numerous and have complicated feedback mechanisms.In combination with additional anthropogenic stressors they can result in a decrease in phytoplankton diversity and dominance of certain species that are more adapted to warmer conditions, such as cyanobacteria and chlorophyta, with cascading effects on zooplankton and the fish stock.The dynamics of tufa formation, a crucial process in the Plitvice Lakes system, would also be influenced by a rise in the lake water temperature.This, in turn, would accelerate the evolution of the lake system's morphology and bring about alterations in its overall structure.
Given the geographical location and climate of Lake Kozjak, along its specific morphology, the findings and conclusions of this study provide valuable insights into the potential response of comparable temperate lakes in mountainous regions.Additionally, the obtained results provide a foundation for future studies in various fields, such as lake biology, geochemistry, sedimentology, and others.

Fig. 1
Fig. 1 Location (a) and a closer view (b) of Lake Kozjak, where the red pin (ϕ = 44.8902°N, λ = 15.6038°E, height of the lake surface 535 m a.s.l.) indicates the location of the maximum lake depth (Source: © Google Earth), and position of Lake Kozjak in the Plitvice Lakes chain (c) (Source: adapted from Klaić et al. 2018)

Fig. 2
Fig. 2 Monthly means of the predicted water temperature vertical profiles for Lake Kozjak in 2019 before (blue) and after (red) model calibration compared to the observed values (black)

Fig. 4
Fig. 4 Yearly average air temperature (a), atmospheric pressure (b), air humidity (c), and wind speed (d) and total yearly precipitation (e) from 2006 to 2100 obtained from the NCC NorESM1-M general circulation model for RCP2.6 and RCP4.5 and the IPSL CM5A-MR general circulation model for RCP8.5 and the REMO2015 regional climate model (r1i1p1 ensemble, v1 downscaling realization) for the Lake Kozjak location; and yearly average shortwave radiation for the same period calculated by SIMO v1.0 (f)

Fig. 5 Fig. 6
Figure5shows the yearly averaged values of the water temperature and the vertical gradient at different depths for the three different scenarios, namely, RCP2.6,RCP4.5, and RCP8.5.While there was an obvious increase in the water temperature under scenario RCP4.5 (Fig.5b1) and even more so under RCP8.5 (Fig.5c1), it was not accompanied by a significant increase in the vertical temperature gradient (Fig.5a2-c2) as the temperature increase occurred at all depths.Furthermore, there is no significant trend of the depth of the maximum vertical gradient indicating no significant trend of the average thermocline depth.As it will be later demonstrated in more detail, significant increase in the maximum thermocline depth is observed under scenario RCP8.5, along with extension of the stratification period in all three analyzed scenarios.However, these findings are not reflected in Fig.5due to the annual averaging.To further examine the stratification evolution, linear regression of the yearly averaged water temperature was calculated for each depth and scenario.Figure6a1-c1show the yearly averaged water temperature in each layer, along with the thermocline temperature.The thermocline temperature is defined as the temperature at the depth of the highest temperature gradient.Since the temperature gradient is calculated at the border between two layers, the thermocline temperature is calculated as the temperature average of those two layers.Since Fig.6a1-c1show yearly averaged temperature values, the thermocline temperature was also calculated as a yearly average using the whole year data and not just data from the periods when the lake was stratified which is generally not the usual approach; so, this should be taken into account if the results are compared to other studies.The thermocline temperature exhibits a significant increase trend of 0.8, 1.82 and 4.68 °C (100 y) −1 under scenarios RCP2.6,RCP4.5 and RCP8.5, respectively.The temperature change rate coefficients were also calculated for each layer and they are shown as a function of the water depth in Fig.6a2-c2.Under the baseline scenario

Fig. 9
Fig. 9 Change rate of monthly values of air temperature (a), total heat flux (b) and heat flux components (cf) calculated for the period from 2006 to 2100 for the RCP2.6,RCP4.5 and RCP8.5 scenarios.Positive heat flux is defined toward the lake surface

Fig. 11
Fig. 11 Water temperature from 2006 to 2010 (a-c) and from 2096 to 2100 (d-f) and Schmidt stability for the same two periods (g-i) for the RCP2.6,RCP4.5 and RCP8.5 scenarios.Dashed lines represent the stratification criterion for the Schmidt stability index for Lake Kozjak of 460 J m. −2

Fig. 12
Fig. 12 Ten-year moving average of the yearly maximum (solid lines) and yearly average (dashed lines) values of the Schmidt stability (a) and the thermocline depth (b) from 2006 to 2100 for the RCP2.6,RCP4.5 and RCP8.5 scenarios