National natural heritage at risk: The Seven Rila Lakes

The area of the Seven Rila Lakes is situated in the North-West part of the Rila Mountain at an altitude between 2100 and 2500 m a.s.l. within the borders of Rila National Park. There are 140 glacial lakes in Rila Mountain but the area of the Seven Rila Lakes is the most visited and famous with its natural beauty and sacred significance. It is a valuable part of the national natural heritage. The popularity of this area grows constantly with the number of visitors from the country and abroad. This process leads to the deterioration of the natural conditions in the site. At the same time, it is not clear what is the effect of climate change on the environment in the lake’s area. There are many factors that contribute to the degradation of ecosystems in the protected area of the Seven Rila Lakes and pose risk to this valuable natural heritage. These factors are natural (mainly climate change-related), anthropogenic (associated with the excessive tourist pressures on the ecosystems in the protected area), and management (stemming from the ongoing conservation policy over the years). This study explores to what extent climate change may put at risk the ecosystems of Seven Rila Lakes. Mean monthly data from The European Center for Medium-Range Weather Forecasts (ECMWF) ERA5-Land reanalysis were used in this study. The resolution of these data is 0.1 x 0.1 (9 x 9 km) and their period is 1981-2020. Reanalysis data include air temperature, precipitation, evaporation, snow depth, etc. Based on data from various sources such as reanalysis data, in situ measurements, and statistical modeling, a scenario, based on current trends in different climatic elements, has been developed in order to project future changes and their impact on lake ecosystems. The results of the modelling of climate change show that in the coming decades an increase in temperature is expected in the high mountain regions of South-western Bulgaria and in particular in the Seven Rila Lakes area. This, combined with the ever-increasing flow of tourists, and high demand for the provided cultural ecosystem services, and insufficient management practices, put at risk the state of the lakes and their capacity to provide the same quality of cultural ecosystem services in the future, which attracts tourists in the area now. Recommendations have been made for the optimization of the management of the protected area in accordance with the observed trends. ABSTRACT


Introduction
The region of the Seven Rila Lakes is a Protected Area under NATURA 2000 in the borders of Natural Park "Rila". According to the Biodiversity Act, there are 30 bird species protected for the purpose of preserving and maintaining their habitats (State Gazette 2008). Only this may be a good enough reason to consider the place a natural heritage, but it is also a territory that has a very high aesthetic value of the landscape, which has given rise to aspirations for spiritual unity with nature among the followers of teacher Dunov, also known as the "White Brotherhood". The lake system is of glacial origin and its catchment is a model for the processes of formation of the glacial type of relief in Southeast Europe through the Pleistocene, which gives it important cognitive significance. The region also contributes to the development of mountain tourism in the country and to economic development of the adjacent municipalities. If to possess an outstanding universal value. However, for Bulgaria, the lakes are undoubtedly subject to national natural heritage, which also meets the accepted definition that "Natural heritage is a geospatial natural element of the socio-ecological system, which provides material and spiritual benefits of timeless importance for previous, present and future generations" (Nikolova et al. 2021c).
A World Heritage property can be entered in on the List of World Heritage in Danger, according to the Article 11 (4) of the Convention, when the condition of the property corresponds to at least one of the criteria in paragraphs 179-180 of the Operational Guidelines. Currently, there are 52 properties included in the List of World Heritage in danger by the World Heritage Committee based on two groups of criteria: Ascertained Danger and Potential Danger. One of the five criteria from the Potential Danger group is the "threatening impacts of climatic, geological or other environmental factors." (World Heritage in Danger, 2021). Out of the List of World Heritage in Danger, there are many threatened sites with cultural or natural significance. We consider that natural heritage sites are at risk when negative changes in their condition lead to the loss of their socioecological value.
In Bulgaria, there is no list of the National Natural Heritage in Danger, but there are examples of such sites. One such example is Lake Srebarna, which has been a biosphere reserve from the World Heritage List since 1983. In the 1970s the lake ecosystem underwent significant changes towards strong eutrophication and anthropogenic speeded-up succession because of the interrupted natural connection with the Danube water due to both natural and man-made factors (Nikolova et al., 2010). There were significant investments in projects to save the lake and as a result, its significance as an Important Bird Area was preserved. Another site from Bulgaria, which is on the World Heritage List, is National Park "Pirin". The huge pressure from the growing tourist infrastructure in the region of the largest ski resort in the country, the town of Bansko, has provoked a series of protests of the public and environmental non-governmental organizations. The other place in Bulgaria that provokes mass protests and active actions to protect it from the ever-growing number of tourists is the Seven Rila Lakes area. Zafirov (2020) provides a detailed analysis of the history of these interactions between civil society and the responsible state institutions regarding the management of the two national parks, "Rila" and "Pirin".
Many unique sites of natural and cultural significance are threatened often by environmental changes all over the World. The International Council on Monuments and Sites (ICOMOS) publishes World Report on Monuments and Sites in Danger (Heritage at Risk) periodically since 2000. In the latest volume (Machat and Ziesemer, 2020), heritage sites at risk are reported from 23 countries, including Bulgaria, in the period 2016-2019. Many of the observed threats in this document relate to climate change. Climate change is a high or very high threat in 33% out of 252 natural World Heritage sites in 2020, while in 2014 they were 15% out of 228 sites listed at that time (Osipova et al., 2020). The main threats to the natural heritage, which are already observed in many sites, are related to the changes of the ecosystem functions, like changes in the timing of biological nutrient cycles, changes in predator-prey, changes in ecosystem services such as pest control, soil stabilization and pollination, an increase of the spread of invasive alien species etc. (Colette Ed., 2007).
Especially sensitive to the pressures of climate change are ecosystems in mountain regions. In the last IPCC report, there is a chapter on the impact of climate change in mountains, which explains the changes which are observed in these environments during recent decades. They include altitude-dependent warming, the general decline of low-altitude snow cover, glaciers and permafrost, change in the precipitation patterns etc. The report also says that "cultural assets, such as snow-and ice-covered peaks in many UNESCO World Heritage sites, and tourism and recreation activities, are expected to be negatively affected by future cryospheric change in many regions" (Hock et al. 2019). The climate change issues in mountain environments were also addressed at the COP 26 Climate Summit. The forum delegates came to an agreement in Glasgow to end deforestation, cut methane emissions, phase-down coal, increase adaptation finance and commit to limit warming to 1.5 degrees by 2030 (Berry 2021).
Climate change could be a threat also to the Seven Rila Lakes in the future, according to the research by a number of authors (Nojarov 2011(Nojarov , 2012a(Nojarov , 2012b(Nojarov , 2012cNikolova 2014;Nojarov et al. 2019Nojarov et al. , 2020Grunewald et al. 2016;Zafirov 2021), as well as according to the regional models of climate change. It is projected with high confidence a decrease of precipitation, lake and river ice, frost, cold spells, permafrost, mean wind speed, and increase of the air temperature, extreme heat, hydrological, agricultural and ecological droughts, and fire weather (IPCC WGI Interactive Atlas, 2021). A rise of the freezing level height and a decrease of snow cover in mountain areas are projected in the future which will impact snow and ice conditions (IPCC AR6 WGI. Regional Fact Sheet -Mountains. 2021). Exploring the local dimensions of these global and regional changes, especially for natural heritage sites in mountainous regions, is a serious challenge. However, they are important for the implementation of effective adaptation models and management of the induced by the climate change risk. Tourism puts more pressure on these sensitive ecosystems. It is especially well expressed during the summer months in the region of the Seven Rila Lakes after the commissioning of the chairlift from "Pionerska" hut to "Rila Lakes" in 2009 (Nikolova et al. 2012(Nikolova et al. , 2012aMitova 2020;Zafirov 2020Zafirov , 2021.
The most visible damage that tourists inflict in this sensitive high-mountain area are pollution, burning of fire, long-lasting camping of many people at the same time, deepening and expansion of the trails. This leads to the destruction of the scarce vegetation, activation of erosion processes and accumulation of deposits in lakes, and contamination of their waters. These adverse effects may be regulated by the application of appropriate management models in a relatively short period of time, while the risk of climate change requires more precise and lengthy regional research. For these reasons, in this study we aim to assess the extent to which the Seven Rila Lakes are threatened by climate change and how these changes will affect the state of the lakes in the future under current conditions and management.

Study area
Rila Mountain belongs to the Rilo-Rhodope Massif. It is the highest mountain on the Balkan Peninsula (Musala peak -2925 m). The region of the Seven Rila Lakes is located in the Northwestern part of Rila Mountain, within the boundaries of the National Park "Rila" (Fig. 1).
This part of the mountain is composed mainly of metamorphic rocks of Pre-Cambrian age, Paleozoic granite, granite-gneiss and Paleogenic sediments. The main quaternary rock materials are result of the exaration, transport and accumulation activity of glaciers and periglacial processes. In this part of the mountain, the slopes go down stepwise to the Valley of the Dzherman River in west and northwest directions. Tectonic features in this part of the mountain have created a dense network of cracks in the root rocks. As a result of the ice age through the Pleistocene, the area has a typical Alpine relief. The seven stepwise circuses in which the lakes are situated were formed through the Würm. The existence of stepped circuses implies repeated migration of the snowline during the ice age. Haramiyata peak is a typical carling developed among the circuses. The movement of the glacier is marked by polished oval sloping rocks, so-called sheep humps, which are well recognised around the lakes Babreka, Bliznaka, and Trilistnika. Тhe moraine ridge supporting the Babreka, Bliznaka and Dolnoto Lakes were formed аs a result of the retreat of the glacier. Sections of trough valley are preserved in the Valley of the Dzherman River downstream the Seven Rila Lakes. An important role in the formation of meso and micro relief forms in the periglacial zone plays the cracked rock base, which contributes to the formation of avalanche spouts and rockfalls shafts. There is a typical 5,5 ha stone glacier of large granite blocks under Haramiyata peak at an altitude of about 2200-2300 m.a.s.l. Dimitrov 2010, Gachev et al. 2017). The processes of solifluction are also widespread. Shallowly developed soils, easily susceptible to erosion, are characteristic of this high mountain area. Peaty marsh soils (Umbric Gleysols), mountain-meadow (Orthic Umbrosols), rankers (Umbric Leptosols) and litosoils (Lithic Leptosols) are common (Mitova 2020). Intense snowmelting, torrential rainfall, and tourist flow during the warm half-year contribute to the entry of large amounts of sediment into the lakes, and these deposits have formed a kind of delta at the infusion of streams from the Okoto and Babreka into the Bliznaka. These conditions characterize the area as extremely sensitive to impacts of both natural and anthropogenic kinds. At the same time, the relief in this part of the mountain has extremely high cognitive and aesthetic value and provides excellent conditions for practicing different types of mountain tourism during the four seasons of the year.
The circus of the Seven Rila Lakes is drained by the River Dzherman. The highest lake is Salzata Lake (2535 m. a. s. l.), followed by Okoto, Babreka, Bliznaka, Trilistnika, Ribnoto and Dolnoto Lake (2095 m. a. s. l.). The lakes are connected by small streams and waterfalls in a lake system. Salzata drains directly into the southern part of Bliznaka, and the streams draining Okoto and Babreka merge into a single stream before its infusion into Bliznaka. The Bliznaka drains into the Trilistnika. The Trilistnika has two drains -one on its eastern side, which flows to Ribnoto Lake, and one at the northwestern end of the lake, which starts the Dzherman River. In the case of low water, through an artificial facility, the current is directed entirely to the Ribnoto Lake. The waters of Ribnoto Lake drain and flow into the Dolnoto Lake, adjacent to the stream that drains it into the Dzherman River.
There are 14 types of habitats with conservation significance in the area (Vasilev et al. 2013). According to EUNIS Classification, the lake system belongs to the habitat class "Oligotrophic to mesotrophic standing waters with vegetation of Littorelletea uniflorae and/or of Isoeto-Nanojuncetea" (EUNIS: C1.1), in category "Vulnerable". "Lakes from this habitat are characterized by high transparency (10-15 m), low temperatures of 10-12°C, blue to dark blue water colour, oxygen composition of 10-14 mg/dm 3 , and neutral acidity pH 6.8-7.2. The Seven Rila Lakes, according to their characteristics and evolution, are in their second or third stage of development, except the lake Okoto which should be in the first stage of development" (Oligotrophic mountain lakes, 2021).
In winter, lakes freeze, and average monthly air temperatures are negative from December to March with a minimum in January. Snowmelt is most intense in May and June when lake levels rise to the maximum extend and become natural reservoirs of snow water. In summer and autumn, the lakes levels decrease depending on temperatures and precipitation, as they do not have significant underground feeding. The results of the measurements in the region of the Seven Rila Lakes in the period 2012-2017 show a good correlation of the air temperature with that of Musala Peak, and with the temperature of the surface lake water (Nojarov, 2015). Over the past decades, from 1991 to 2020, the average annual air temperature in the country was 11.8°C. The norm for the annual amount of precipitation above 2000 m altitude is 850 mm (NIMH, 2021). There are statistically significant trends (per decade) of 0.5 °C for the mean annual temperatures and of 0.7 °C for the mean monthly temperature (May-September) during the period 1994 -2017 at the Musala peak. The trend of the annual precipitation amount (mm) in the same period is -5.8 but it is not statistically significant (Nojarov et al. 2019). Rising temperatures and decreasing precipitation would lead to a drop in the levels and a corresponding reduction of the surface areas of the Seven Rila Lakes. Remote sensing surveys have identified such changes, and rise of the upper limit of the forest in this part of the mountain (Gikov and Dimitrov 2010, Nikolova et al. 2012, Zafirov 2020).
There are comparatively fewer studies on the ecosystem services and benefits provided by the Seven Rila Lakes area for the development of tourism and recreation, as well as those for assessing the balance between the recreational potential of this territory and the impact of tourism on it. In 2014, the first assessment and mapping of the ecosystem services in the area of the Seven Rila Lakes was carried out by Nedkov et al. (2014). The results of this study show that the ecosystems in the area of the Seven Rila Lakes provide provisional, regulating and cultural services. "The water bodies have very high cultural and provisional capacity, and although their regulation capacity score is lower, they should be identified as the main object of protection" (Nedkov et al, 2014).
The studies on the adverse effect that tourism has on the condition of the protected area have started after the commissioning of the chairlift from "Pionerska" hut to hut "Rila Lakes" in 2009 Nikolova et al. (2012, 2012a) Mitova (2020, Zafirov (2020Zafirov ( , 2021, etc. There are many media and NGO publications on the impact of tourism pressure on the lakes' system, as well as some protests in defense of the lakes which were analyzed in details in Zafirov (2020).

Data and methods
The data for the Seven Rila Lakes were collected through expedition measurements carried out during the warm half of the year in the period 2012 -2020. The total number of expeditions, respectively measurements, is 21. Water levels and temperatures of the seven lakes were parameters being monitored. At the beginning of the measurements (28.06.2012), a marker was placed on rocks in each of the lakes, which marks level 0. All subsequent measurements of water levels are against these markers. The water temperature was measured at a depth of 20-30 cm near the shore. Measurements of water levels were carried out on all expedition dates (21) except on one date for each of the lakes Dolno, Ribno, Okoto and Salzata, due to bad weather conditions. There are three missing measurements of water temperature (respectively a total of 18 cases) for the lakes Dolno, Ribno, Trilistnika, Bliznaka and Salzata and four missing measurements (respectively a total of 17 cases) for the lakes Babreka and Okoto. Additionally, the data of the first (28.06.2012) measurement of water temperature of the lakes Okoto and Salzata were taken out of the study, because at that moment their ice cover was not yet completely melted, and therefore the presence of ice in the water had a significant impact by lowering the water temperature.
Hourly and monthly averaged data from The European Center for Medium-Range Weather Forecasts (ECMWF), ERA5-Land reanalysis (Copernicus Climate Change Service (C3S) 2017) were used for the studied region. The resolution of these data is 0.1 x 0.1° (9 x 9 km) and they cover the period 1981 -2020. The gridcell, which completely covers the area of the seven Rila lakes, has coordinates 42.1-42.2°N and 23.3-23.4°E. For the purposes of the present study, the data for the following elements were used -air temperature at 2 m above the surface, precipitation, evaporation, U-wind at 10 m above the surface, V-wind at 10 m above the surface, snow depth water equivalent, runoff and relative humidity (RH). The water balance of the studied grid-cell was also calculated according to the formula: where B is water balance, P is precipitation, SD is snow depth water equivalent, E is evaporation (having always negative values in ERA5-Land reanalysis) and R is runoff.
The warm half-year in high mountains of Bulgaria (Rila and Pirin) generally includes the period May-October (Nojarov et al. 2019). During this period the lakes in the study area are partially or entirely free of ice. Because water levels of the lakes on a certain date are a function of the entire previous period, in studying and modeling of water levels the data for the different climatic elements should be averaged for the entire period from the beginning of May to the date of the respective measurement. If the measurement date is at the beginning of the month, average monthly values of the previous months (starting with May) are averaged and used in further calculations. If the date of measurement is at the end of the month, both average monthly values of the previous months (starting with May) and average value of the month of measurement are averaged and used in further calculations. If the measurement date is in the middle of the month, average monthly values from the previous months (starting with May) and average daily values (based on the hourly data from ERA5-Land reanalysis) up to the specific day of measurement are averaged and used in further calculations. Since the temperature of the water is relatively quickly responsive to environmental conditions, in its studying and modeling, average monthly values of the different climatic elements for the month, in which the respective measurement was carried out, were used.
This study uses mainly statistical methods (Wilks 2006). The level of statistical significance for all calculations is p < 0.05. The trend analysis was done by means of linear regression. The linear regression represents the relationship between the independent variable (time) and the dependent variable (air temperature, precipitation, evaporation, U-wind, V-wind, etc.) through a linear equation based on the observed values. The trends of different climatic elements were calculated for the warm (May -October) half of the year for the period 1981 -2020. Spearman's rankorder correlation was used in order to establish the relationships between the various climatic elements, and lake levels and water temperatures. The advantage of this correlation is that it does not depend on the type of statistical distribution of the studied variables. Statistical distributions of the variables were also checked and the results showed that they are close to normal distribution. This makes it possible to use multiple linear regression models to show the complex influence of different climatic elements on lake levels and water temperatures (following, generally, the methodology described in Nojarov et al. (2020) and Nojarov et al. (2021a). The climatic elements serve as potential predictor variables in the created statistical models, and lake levels and water temperatures are predicted variables. The statistical models are created separately for lake levels and water temperatures of the lakes. The method which was used in multiple linear regression is forward stepwise method. It adds predictors one at a time as long as these additions are worthy (the inclusion should give the most statistically significant improvement of the fit). After a predictor has been added, all predictors in the current model are checked to see if any of them should be removed. Then the process continues until a stopping criterion is met. The criterion that was used as to which predictor (variable) would eventually be included in the statistical model was minimum Corrected Akaike Information Criterion (AICC). It prevents well an overfitting of the created multiple linear regression model in case of too many potential predictors. These multiple linear regression models are then used to calculate the projected values of lake levels and water temperatures for the next 50 years (2021-2070). These calculations use as an input the projected linear regression values of the predictors that have statistically significant trends for the period 1981-2020. Those variables, which do not have a statistically significant trend for the period 1981-2020, are included in these calculations with their average values for the period 1981-2020 (i.e., no projected change in their values).

Results
The trends during the warm half of the year (May -October) of the different studied climatic elements in the region of the seven Rila lakes for the period 1981 -2020 are shown in Table 2. A statistically significant positive trend is observed in air temperature. To a large extent, evaporation is related to air temperature, which explains the observed statistically significant negative trend. It should be borne in mind that evaporation is measured in negative values and a negative trend actually means an increase in evaporation. There is a statistically significant negative trend in U-wind as well.
Considering that the positive values of this element represent a wind with a westerly direction and vice versa, the observed trend means that in the area of the seven Rila lakes in recent decades the transport of air masses from the east is increasing during the warm half of the year. A statistically significant negative trend is observed in snow depth water equivalent. In general, snow cover exists at the beginning of the warm half of the year, usually in May and June, and it represents the conditions in the cold half of the year. It also increases the amount of water in the lakes due to its melting at the beginning of the warm half of the year. The observed negative trend in this element is mainly related to the rising air temperatures, which leads to an earlier melting of the snow cover accumulated during the cold half of the year. This negative trend also affects the water balance of the studied area, which also decreases with a statistically significant trend. Snow depth water equivalent is one of the important positive terms of the water balance equation. The other important positive term is precipitation, which, however, has not changed significantly in the period 1981-2020.
The coefficients of the Spearman correlation between water levels of the seven Rila lakes and the different climatic elements are shown in Table 3. The period here is 2012-2020. The values of the climatic elements are calculated as explained in the data and methods National natural heritage at risk: The Seven Rila Lakes section. The results show that air temperature is negatively correlated with water levels of the lakes, but the coefficients are statistically significant only for the two highest lakes -Okoto and Salzata. In general, both air temperature and evaporation relationships with water levels are mediated by precipitation. Usually in meteorological conditions with precipitation, air temperature and evaporation are lower, and at the same time, the water levels in the lakes rise. Thus, it can be summarized that the coefficients for air temperature and evaporation in Table 3 represent an indirect relationship. The relationship between precipitation and water levels of the lakes is direct, which is evidenced by the fact that all coefficients in the relevant column of Table 3 are positive and statistically significant. The correlation coefficients between U-wind and V-wind on one hand and water levels of the lakes on other hand also have positive values and are statistically significant for most of the lakes (except for Okoto). This relationship is also largely due to precipitation, because the transport of air masses from south-west during the warm half of the year is associated with more precipitation, and the transport of air masses from north-east is associated with less precipitation. The runoff is in a statistically significant positive correlation with water levels because higher levels are associated with more runoff from the studied area and vice versa. The snow depth water equivalent is in a positive but statistically insignificant (except for Okoto) relationship with water levels of the lakes. It is a positive term in the water balance equation and should higher water levels of the lakes as it melts, but in recent years as a result of warming snow cover in the warm half of the year is less and less, and has become an insignificant factor for water levels. The relative humidity has statistically significant positive correlation coefficients with water levels of the lakes. This relationship is also to some extent mediated by precipitation, because in meteorological conditions with clouds and precipitation, the relative humidity is higher and vice versa. The water balance of the territory is positively correlated with water levels, but plays more significant role in the four higher lakes. This can be explained to a large extent by the anthropogenic intervention in the lower lakes, the water of which is used for various utility needs, which affects their water levels.
The coefficients of Spearman correlation between the water temperature of the seven Rila lakes and the different climatic elements are shown in Table 4. The period here is also 2012 -2020. The values of the climatic elements are calculated as explained in the data and methods section. The climatic elements evaporation, runoff, and water balance are not considered here, because they do not affect water temperature. As could be expected, the best relationship with water temperature has air temperature. The correlation coefficients are high, and statistically significant and positive for all lakes. Precipitation and U-wind do not show any significant relationship with water temperature. V-wind has a statistically significant negative correlation with water temperature. This means that when the wind is from the north, the water temperature is higher and vice versa. Usually, the transport of air masses from the north during the warm half of the year is associated with more stable, anticyclonic and cloudless weather in the Bulgarian mountains, which leads to higher air temperatures, and hence to higher water temperatures. The snow depth water equivalent is significantly negatively correlated with water temperature, because  melting of the snow leads to an inflow of colder water into the lakes and accordingly lowers their water temperature. This means that the more snow there is to melt in the warm half of the year, the lower the water temperature of the lakes is and vice versa. Relative humidity has a statistically significant negative relationship with water temperature. This is again related to the specific meteorological conditions, when in cloudy, rainy, and with high relative humidity weather the air temperature is usually low and this leads to lower water temperatures. Thus, this relationship is not direct but is mediated by the air temperature. It can be concluded that when the relationships of the different climatic elements with water levels and temperatures of the lakes are considered separately, the main factor, which has direct influence on water levels is precipitation, and the main factor, which has direct influence on water temperatures is air temperature. Snow depth water equivalent in the warm half of the year is a secondary factor that has a direct influence on water levels and temperatures of the lakes Table 5 shows the results of the multiple linear regression models, the dependent variables being water levels and water temperatures of the seven Rila lakes (first column of the table). The table also includes the values of adjusted R 2 , which show how well the statistical model describes variations of the dependent variable, as well as the predictors in the respective model with the sign of their relationship with the dependent variable (plus or minus sign) and the percentage of variation of the dependent variable that they describe. These MLR models are based on data from measurements in lakes in the period 2012-2020. The results in Table 5, regarding the values of adjusted R 2, show that water levels of the lakes are modeled worse than the temperature of the water. This is mainly due to two reasons. The first one is that the lower lakes (Dolno, Ribno, and Trilistnika) are subject of a significant anthropogenic intervention -the use of lake waters for utility needs, as well as for electricity production. Thus, the influence of natural climatic factors on the water levels of the lakes decreases. This is also reflected in the lower values of adjusted R 2 . The second reason, which is valid for all seven studied lakes, is the way the measurements are taken, on an expeditionary basis on certain dates during the warm half of the year. In general, water levels react relatively quickly to precipitation. This means that if the measurement is made immediately after a precipitation spell, the water level will be higher. Accordingly, if the measurement is made long after the last precipitation, the water level will be lower. The idea of this study is to reveal the impact of climatic elements on water levels for a longer term, i.e. for the entire warm half of the respective year. That is why the data of the predictors included in the MLR model are calculated on this basis. Short-term fluctuations of water levels are not outputted by the models, however, they affect the measurements, which explains the low values of adjusted R 2 . Thus, the results obtained here, as well as projected values regarding water levels of the seven Rila lakes should be understood as water levels at the end of the warm half of the year (September, October) under conditions of average precipitation in the entire warm half of the year.
Comparisons between the values of the water levels of the seven lakes modeled by the statistical models with measured (observed) values for the period 2012-2020 are shown in Figure 2. It can be seen that, in general, water levels are better modeled after 2014, which is due to the fact that after this year the frequency of measurements decreases and the possibility of a certain measurement to be made immediately after precipitation event also decreases, but does not disappear completely. Such is the case of the last measurement in October 2020, when there is a significant discrepancy between the modeled and measured values because of a precipitation spell before the date of the measurement (except for Okoto, which is the deepest and has the largest water volume lake and the change in its water level relative to given precipitation is the smallest, as can be seen in Table 5, where the value of adjusted R 2 is the highest). The predictors that explain the variations in water levels of the lower five lakes are precipitation (Dolno, Ribno) and U-wind (Trilistnika, Bliznaka, and Babreka). Higher precipitation and U-wind values lead to higher water levels and vice versa. The positive values of U-wind represent wind from the west, which during the warm half of the year in Bulgaria is usually associated with cyclonic conditions and, respectively, precipitation. There are other predictors that determine water levels of lakes Okoto and Salzata. For Okoto, these are runoff and V-wind, which are in a positive relationship with the water level. It is clear that higher water levels are associated with more runoff and vice versa, while for V-wind a positive relationship means that wind from the south is associated with higher water levels and vice versa. In general, the wind from the south is associated with cyclonic Table 5. The results of the multiple linear regression models for water levels and water temperatures of the seven Rila lakes (first column). The table also includes the values of adjusted R 2 and the predictors in the respective model with the sign of their relationship with the dependent variable (plus or minus sign) and the percentage of variation of the dependent variable that they describe. These multiple linear regression models are based on data for the period 2012-2020. RH stands for relative humidity. RH+12.1% V-wind-8% V-wind-9.6% V-wind-9.9% V-wind-9.4%  conditions over Bulgaria and, accordingly, with more precipitation during the warm half of the year. For Salzata, the predictors are water balance and relative humidity, which are positively correlated with the water level. The relative humidity, as mentioned earlier, is in directly proportional relationship with precipitation in mountainous territories of Bulgaria (Nojarov 2012c), which explains the sign of the relationship in the statistical model. The water temperature of the seven Rila lakes is modeled very well by MLR models. This is evident both from the values of adjusted R 2 in Table 5 and from the graphs in Figure 3, which show the modeled and observed values of water temperature. The main predictor for all lakes is air temperature, which is in directly proportional relationship with water temperature. These results are in agreement with the results in Table 4. The temperature of the surface layer of the water in the lakes reacts relatively quickly to changes in air temperature. Another predictor that is frequently present is V-wind, which is in an inversely proportional relationship with water temperature. This means that wind from the north (negative V-wind values) leads to higher water temperatures and vice versa. In general, wind from the north in the Bulgarian mountains during the warm half of the year is associated with anticyclonic conditions, sunny weather, and higher air temperatures.
The described above statistical models were used to project the future (2021-2070) tendencies in water levels and temperatures of the seven Rila lakes. The results of the projections of water levels of the lakes for the period 2021-2070 are shown in Figure 4. Once again it should be noted that these are the water levels at the end of the warm half of the year (September, October) with an average precipitation for the entire warm half of the year. The figure shows that for four of the lakes (Trilistnika, Bliznaka, Babreka, and Salzata) a decrease in water levels is expected, which is in the range of 8-10 cm by 2070 for three of the lakes and of about 1 cm by 2070 for Salzata. For the other three lakes (Dolno, Ribno, and Okoto) no change in their water levels are expected. In general, based on these results and the uncertainty of the modeling, it can be concluded that the water levels of the seven Rila lakes over the next 50 years will remain unchanged or there will be a slight decrease.
The results of the projections of water temperature of the seven Rila lakes for the period 2021-2070 derived from the statistical models are shown in Figure 5. It reveals that the water temperature  of all lakes will increase by 2070. This increase will be about 1°С for Okoto and Salzata, about 1.5°С for Trilistnika, Bliznaka and Babreka and about 2°С for Dolno and Ribno. The increase in water temperatures of the lakes is mainly due to the expected increase in air temperatures in the mountains of Bulgaria over the next 50 years.

Discussion
The results in the present study show that the trends in the atmosphere over the seven Rila lakes do not differ from those existent in the wider territory of Bulgaria or Southeastern Europe. A statistically significant positive trend is observed in air temperature, which is in line with the results from other studies on the mountains of Bulgaria (Nojarov 2011, Nojarov 2012a, Nojarov 2012b, Nojarov 2015, Nojarov et al. 2019. Also, as a result of global warming, there is a change in atmospheric circulation during the warm half of the year as the transport of air masses from the north and east increases, which corresponds to the results obtained by Nojarov (2021b). A similar process is observed also in Pirin Mountains (Nojarov et al. 2019).
Projected values of the lake's water temperature increase most in the lakes situated between 2095 and 2216 m a.s.l. and much less for the lakes situated between 2440 and 2535 m a.s.l. This reflects not only the direct connection with air temperature, but also relates to the lake's morphometry. Such changes in the water temperature of lakes will undoubtedly lead to changes in their ecological status. It is the subject of another expert's report to indicate how lake ecosystem would respond to these changes, but they are sufficiently substantial to take appropriate adaptation and mitigation measures.
The climate change impact on the lakes could be accelerated significantly under the current tourism pressure in the studied area. After the chairlift from "Pionerska" hut to "Rila Lakes" hut started its operation in 2009, the number of visitors increased from 15 000 in 2008 to 60 000 in 2009 and to over 100 000 in recent years (Nikolova et al. 2012a, Mitova 2020, Zafirov 2020, Zafirov 2021. The chairlift has a length of 2136 m and a capacity to carry 960 people per hour, according to data from National Park "Rila" (2012). The better accessibility to the lakes leads to a drastic increase of the anthropogenic pressure. The exceptional natural beauty of this place affects tourists very strongly and emotionally and has a lasting impression on them. According to a survey carried out in 2012, 33.5% of 63 respondents identified the area of the Seven Rila Lakes as "the most beautiful place" in Bulgaria and another 26.4% as "extremely beautiful". All respondents recommended to their friends to visit the lakes (Nikolova et al. 2012a). These results confirm the high capacity of the area to provide cultural ecosystem services for tourism and recreation estimated by Nedkov et al (2014). At the same time, they show how fast increase the demand of the provided cultural services. Mitova (2020) investigate the trails and campsites impacts on the Seven Rila Lakes. The results show that the visitor flow has created a trails network with an average density of 75.4 m/ ha in the studied area of 457.1 ha. Two campsites are located on 1.4% of this territory during the summer tourist season. About 5% of the length of the trails exerts strong and extremely strong pressure on soils and geo-heritage. About 14% of the length of the trails has a strong and extremely damaging impact on natural habitats and biodiversity. The study found that the catchment area of Trilistnika Lake is the most threatened by erosion processes caused by tourism and recreation activities.
Data from remote sensing surveys show that the eutrophication of lakes increases progressively during the warm half-year. More than 50% of the water mirror of Ribnoto Lake was occupied by vegetation in September 2010, while in the images from October 1977, the water mirror is transparent and free of vegetation (Nikolova et al. 2012).
An increase in eutrophication was also found in the lakes Ribnoto during the summer months of the period 2017-2020, (Zafirov, 2021). The increase of water temperature and the nitrogen loading from the eroded soils intensify the process of eutrophication and lake deterioration.
The management of the National Park "Rila" is in accordance with the Protected Areas Act, the Biodiversity Act, the Management Plans of the Rila National Park, and other regulations. Among the recommended conservation measures for the Seven Rila Lakes given in the National Action Plan for Conservation of Wetlands of High Significance in Bulgaria (2013) is "Regulation of the tourist flow and all tourist-related services" and "Specialized monitoring of the succession processes and, if necessary, artificial limitation of competitive species". The Ministry of Environment and Water and the administration of National Park Rila implemented measures for the reduction of the anthropogenic pressure in the lake's area. The Minister of Environment and Water, by order from March 2021 established a protected area "Rila Buffer" (State Gazette 2020). Its border goes down to the starting station of the lift. The access to the Seven Rila Lakes with high-passable cars for tourism purposes was banned and the project to build hiking trails in the lakes area was renewed in 2020 in order to reduce the dispersed movement of tourists in the lakes area. Each of the last two measures was met with sharp critical reactions from citizens and environmentalists because of different views on how to solve the problems. In 2017 and 2018, the vegetation from Ribnoto Lake was buried by mechanical means. This act, expectedly, does not lead to the desired result, and after about a year the vegetation again filled the lake. According to data from remote surveys by Zafirov (2021), in 2019 and 2020 the values of the NDVI index for the lake even exceed those from before the intervention. The development of each ecosystem is accompanied by a natural process of succession, which can be significantly accelerated by changes in the environment.
Unfortunately, the protected areas in Bulgaria are still rarely seen as natural heritage sites. This would make a significant difference to both their preservation and promotion and their capitalisation as tourist sites. Studies of Nedkov et al. (2014), Nedkov et al. (2018Nedkov et al. ( , 2021Nedkov et al. ( , 2021aNedkov et al. ( 2021b, Zhiyanski et al. (2021), Hristova and Stoycheva (2021), Neugarten et al (2018) provide a good methodological basis for the development of sustainable models for tourism development. The practical application of the ecosystem approach for assessing the capacity of ecosystems to provide ecosystem services for recreation and tourism is well developed both for the territory of the country (Ihtimanski et al. 2020, Prodanova 2021) and for specific natural heritage sites in Bulgaria (Silvestriev et al. 2021, Semerdzhieva and Borissova, 2021).

Conclusion
The main climatic factors that exert influence on water levels of the Seven Rila Lakes are precipitation and snow depth water equivalent. The factors that exert influence on the water temperatures of these lakes are air temperature and snow depth water equivalent. The other studied climatic elements have no direct influence. Their influence is realized either through precipitation or air temperature. Over the last 40 years, there have been trends of a significant increase in air temperature and evaporation during the warm half of the year in the studied area. There is a significant change in atmospheric circulation because the transport of air masses has changed from predominantly westerly to predominantly easterly. The snow cover, which remains unmelted at the beginning of the warm half of the year, is decreasing more and more. This also affects the water balance of the territory, which also has decreased significantly in the period 1981-2020. These trends, combined with the correlations revealed in the study, lead to the conclusion that in the next 50 years water levels of the seven Rila lakes are expected to remain unchanged or to decrease slightly, and their water temperatures in the warm half of the year are expected to rise significantly by about 1-2°C.
The combined impact of these two risk factors, tourist flow, and climate change, accelerates the risk for the current ecological state of the Seven Rila Lakes. The management of the protected area is faced with both the problem caused by tourist pressure and the problem emerging from climate change. More targeted actions are necessary to limit the number of tourists in line with the capacity of this natural heritage site to provide services for recreation and tourism. The management of protected areas also as sites of national natural heritage would significantly expand the tools for their sustainable use and risk mitigation. Climate change is a long-term challenge which needs effective implementation of the National Climate Change Adaptation Strategy and Action Plan in the protected area's management plans.

Funding program
This research was funded by the BG05M2OP001-1.001-0001 Project "Creation and development of "Heritage BG" Centre of Excellence", Operational Programme Science and Education for Smart Growth, Priority axis 1, Procedure BG05M2OP001-1.001, Component 4 "New technologies in creative and recreation industries".