The Combined Impact of Hydropower Plants and Climate Change on River Runoff and Fish Habitats in Lowland Watersheds

: Aquatic ecosystems are particularly vulnerable to anthropogenic activity and climate change. The changes in ﬂow regimes in Lithuanian lowland rivers due to the operation of hydropower plants (HPPs) and the impact of altered ﬂow on some ﬁsh species have already been studied. The impact of climate change on future natural river runoff and the structure of ﬁsh assemblages was also investigated. However, it is still unknown how the combined effect of climate change and ﬂow regulation related to hydropower generation may affect ﬁsh assemblages in the downstream river reaches below the Lithuanian HPPs. In this study, the physical habitat modelling system MesoHABSIM was used to simulate spatial and temporal changes in aquatic habitats availability for different ﬁsh species under the inﬂuence of HPP at different climate change scenarios. Changes in the available habitat were assessed for common ﬁsh species in four HPP-affected rivers representing different hydrological regions of Lithuania. The modelling results showed that the operation of HPP under climate change conditions in most rivers could be beneﬁcial for small benthic ﬁsh species such as gudgeon Gobio gobio and stone loach Barbatula barbatula . Meanwhile, for larger ﬁsh species (e.g., chub Squalius cephalus and vimba Vimba vimba ) the alteration in the temporal availability of suitable habitat was relatively higher. The simulation results show that depending on the climate change scenario and the river region, the impact of the HPP on temporal habitat availability for ﬁsh will become even more signiﬁcant. Already in the near future, this impact should be most noticeable in the HPP-affected Lithuanian rivers of the Central region, which discharge is dominated by surface runoff. In the groundwater-fed rivers of the South-eastern region, a signiﬁcant increase in the impact of HPP on temporal availability of suitable habitats for certain ﬁsh species is more likely in the far future, while in the near future, the possibility of a signiﬁcant impact is higher only at the RCP2.6 climate change scenario. The exception is the rivers in the Western region (represented by the Bartuva River in this study), where small changes in the impact of HPP on ﬁsh habitat are predicted. Only the temporary availability of habitat suitable for a schneider may deviate depending on the climate change scenario.


Introduction
In the modern world, humankind is changing the natural environment at an unprecedented rate. At the beginning of this century, the Millennium Ecosystem Assessment [1] stated that over the past 50 years, humankind has changed Earth's ecosystems more rapidly and extensively than in any comparable period of time in human history. As a result of this ongoing change, more than 1 million plant and animal species are threatened with extinction [2]. The biodiversity of freshwater ecosystems is declining much more than that of the most affected terrestrial ecosystems [3]. Aquatic ecosystems are particularly vulnerable to climate change because (i) many species within these fragmented habitats have limited abilities to disperse as the environment changes; (ii) water temperature and habitat availability are climate-dependent; and (iii) many systems are already exposed to numerous anthropogenic stressors [4].
Although it is common to mention rising water temperatures when it comes to the effects of climate on river ecosystems, changes in flow variability are no less important but more difficult to predict. Projections suggest that climate change would lead to significant changes in the hydrological behavior of rivers in many countries [5][6][7][8]. Climate change has been estimated to significantly alter seasonal flow regimes in 90% of the world's land area [9]. Climate-induced changes in flow regimes cause additional stress to river systems that are already severely affected by anthropogenic activities such as hydropower production and dams [10,11]. The EU Member States point out that hydropower and dams are major contributors to the deterioration of the aquatic environment [12,13]. Forty-eight percent of rivers (expressed as river volume) globally are moderately to severely impacted by either flow regulation, fragmentation, or both [14], yet the hope for the sustainable development of dams remains [15].
Flow variability is considered a major factor in ecological integrity that determines the size, shape, structure, and dynamics of river ecosystems [14,16]. To represent biologically relevant streamflow features, Olden and Poff [17], found as many as 171 hydrological indices used to describe flow regimes in scientific literature adequately. A natural flow paradigm states that the full range of natural intra-and interannual variation of hydrological regimes, and associated characteristics of timing, duration, frequency, and rate of change, are critical in sustaining the full native biodiversity and integrity of aquatic ecosystems [18]. Changes (of any origin) in any of these characteristics can affect key processes in aquatic ecosystems, reduce their resilience and, in the long run, lead to irreversible alterations [19][20][21][22][23]. A number of recent studies have indicated the impact of modified flows on various species of aquatic ecosystems: from primary producers diatoms and macrophytes [24], butterflies, grasshoppers, vascular plants [25], to birds [26]. However, much of the literature pays particular attention to the multiple impacts of altered flow characteristics on fish communities [27][28][29][30][31][32][33][34][35], because fish are the most affected by hydrological disturbances among other groups of aquatic organisms [36,37].
The changes in flow characteristics in Lithuanian lowland rivers due to hydropower generation [38] and the impact of altered flow on some fish species at different annual runoff [39] have already been studied. The impact of climate change on future river runoff [40,41], its extremes [42], and the structure of fish assemblages [43] were also assessed. It has been established that an increase in water temperature at the end of the century will become significant pressure, which will lead to a decrease in the number of stenothermal fish, while eurythermal fish will not be affected [43]. However, it is still unknown how the combined effect of climate change and flow regulation for hydropower generation may affect fish assemblages downstream of hydropower plants (HPPs) in Lithuania. The study seeks to project changes in hydroclimatic and hydromorphological conditions downstream of river dams and to explore how these habitat modifications may affect lowland river fish communities.

Study Area and Data
Lithuania is located in the south-eastern part of the Baltic Sea region. Around 22.2 thousand rivers and small streams lay down within the country territory of 65.3 thousand km 2 . However, as many as~80% of the entire river network consists of streams shorter than 3 km. The rivers of Lithuania distinguish by their main three hydrological responses to local physico-geographical features; therefore, Lithuania is divided into three main hydrological regions-Western, Central, and South-eastern. The Western hydrological region is described as the wettest one because the western slopes of Samogitian Uplands are a primary barrier for the air masses and they collect the highest amount of annual precipitation (750-900 mm), accordingly all rivers of this region are mostly fed by the rainfall. In the Central hydrological region, soils of heavy properties dominate. This feature determines weak infiltration and supplies part from the groundwater. Moreover, the favorable conditions for the formation of snow cover during the cold period are inherent. Due to these reasons, surface feeding, such as snowmelt and rainfall, dominates in the rivers of the Central hydrological region and determines high variability in the flow during the year. The South-eastern hydrological region falls within the Baltic Uplands, where sandy soil properties are common. Such nature of the physical environment significantly affects the behavior of river hydrology by the effective absorption of surface runoff surplus to groundwater during the spring and gradually supply in the summer low-flow period. Due to a general slight elevation of the landscape, the vast majority of Lithuanian rivers can be classified as lowland rivers. River damming produced several hydromorphological impacts on Lithuanian rivers. At least 20 gauged and 36 ungauged rivers are affected by 102 small hydropower plants (Figure 1). These numbers highlight the necessity to understand the impact of HPPs on the hydrological behavior of rivers of the different hydrological regions. Therefore, four river catchments (Verknė, Širvinta, Šešupė, and Bartuva) were chosen to cover each hydrological region. It will help to get a broader picture of the effects of HPPs on the fish habitats and to shift the estimated patterns to ungauged river catchments. the Central hydrological region and determines high variability in the flow during the year. The South-eastern hydrological region falls within the Baltic Uplands, where sandy soil properties are common. Such nature of the physical environment significantly affects the behavior of river hydrology by the effective absorption of surface runoff surplus to groundwater during the spring and gradually supply in the summer low-flow period. Due to a general slight elevation of the landscape, the vast majority of Lithuanian rivers can be classified as lowland rivers. River damming produced several hydromorphological impacts on Lithuanian rivers. At least 20 gauged and 36 ungauged rivers are affected by 102 small hydropower plants ( Figure 1). These numbers highlight the necessity to understand the impact of HPPs on the hydrological behavior of rivers of the different hydrological regions. Therefore, four river catchments (Verknė, Širvinta, Šešupė, and Bartuva) were chosen to cover each hydrological region. It will help to get a broader picture of the effects of HPPs on the fish habitats and to shift the estimated patterns to ungauged river catchments. Hydrometeorological information such as daily discharge data (Q, m 3 /s), daily air temperature (T, °C), and daily precipitation amount (P, mm) was taken from hydrological and meteorological yearbooks of the Lithuanian Hydrometeorological Service. The available discharge data of the following four water gauging stations (WGS): Verknė-Verbyliškis (1952-2018), Širvinta-Liukonys (1972-2018), Šešupė-Kudirkos Naumiestis , and Bartuva-Skuodas (1957 was used for the calculations of different hydrological characteristics as well as calibration and validation of hydrolog- Hydrometeorological information such as daily discharge data (Q, m 3 /s), daily air temperature (T, • C), and daily precipitation amount (P, mm) was taken from hydrological and meteorological yearbooks of the Lithuanian Hydrometeorological Service. The available discharge data of the following four water gauging stations (WGS): Verknė-Verbyliškis (1952-2018), Širvinta-Liukonys (1972-2018), Šešupė-Kudirkos Naumiestis , and Bartuva-Skuodas (1957 was used for the calculations of different hydrological characteristics as well as calibration and validation of hydrological models of the selected rivers. Additionally, the creation of hydrological models requires meteorological information, therefore daily air temperature and precipitation amount of nine meteorological stations (Klaipėda, Telšiai, Tauragė, Lazdijai, Kaunas, Vilnius, Varėna, Utena, and Ukmergė) for the period 1986-2005 was used. Characteristics of selected hydropower plants (HPPs) were taken from the Rules of Reservoir Exploitation and Maintenance (Table 1). Output data of daily air temperature and precipitation amount of RCA4 regional climate model (RCM) based on three driving global climate models (EC-EARTH, HadGEM2-ES, and MPI-ESM-LR) were used for the future climate projections. The projections were simulated according to three RCP scenarios (optimistic-RCP2.6, realistic-RCP4.5, and pessimistic-RCP8.5) for the near (2021-2040) and far (2081-2100) future, as well as for the reference period (1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005). The daily output of meteorological data for the selected combinations of a regional climate model was extracted for 11 × 11 km grid cells from the EURO-CORDEX database (www.euro-cordex.net, accessed on 5 December 2021).
For the meso-habitat modeling, the State river monitoring database was used to select the target fish community of the considered river reaches. The State river monitoring covers data on the physical, chemical, and hydrological characteristics of river sites, as well as fish data collected by standardized single-pass electric fishing in mid-July-September on river sections with a minimum length of at least 10 times the wetted width (but not less than 50 m). Only river sites in natural conditions (from good to high ecological status according to the European Water Framework Directive) with a catchment area >100 km 2 and sampled by wading were selected from the database (in total, 245 sites). Further, fish species were selected, the frequency of occurrence of which in the species-specific types of rivers is >40%, and the relative abundance is >5%. For anadromous species, only the frequency of occurrence criterion was applied. In total, 12 fish species met the criteria and were selected as adequate target fish community: bleak Alburnus alburnus (standard length (SL) 7-13 cm), schneider Alburnoides bipunctatus (SL 9-13 cm), stone loach Barbatula barbatula (SL 10-12 cm), bullhead Cottus gobio (SL 8-10 cm), gudgeon Gobio gobio (SL 10-15 cm), dace Leuciscus leiciscus (SL 20-25 cm), minnow Phoxinus phoxinus (SL 6-9 cm), roach Rutilus rutilus (SL 15-25 cm), juvenile salmon Salmo salar (SL < 15 cm), juvenile trout Salmo trutta (SL < 15 cm), chub Squalius cephalus (SL 25-40 cm) and vimba Vimba vimba (SL 30-50 cm). The listed standard lengths of fish species of Lithuanian inland waters are described according to [44].

Methods
The applied methodological scheme involved different approaches, which together allowed us to evaluate the potential changes in indices of fish habitats under the influence of hydropower plants and climate change ( Figure 2). The first level included the fundamental knowledge on regional differences to select reasonable pilot rivers and to evaluate their potential changes in different hydrological environments. The second level covered the initial data required for the implementation of such kind of research. It consisted of the primary data available from existing databases and datasets, as well as the data that was Water 2021, 13, 3508 5 of 21 necessary to be collected during the field surveys. The third level was composed of the analysis of different kinds of data, i.e., • an analysis of hydrological and hydromorphological data was collected during the field surveys; • development of conditional habitat suitability criteria according to the data of field surveys and fish monitoring; • evaluation of reference and altered (under influence of HPP) conditions according to historical observations and method of analogy; • simulation of the projections of pilot rivers using created HBV hydrological models that required the observed hydrometeorological data for calibration and validation, and output of ensemble of regional climate models according to different RCP climate scenarios for the discharge projections.
the initial data required for the implementation of such kind of research. It consisted of the primary data available from existing databases and datasets, as well as the data that was necessary to be collected during the field surveys. The third level was composed of the analysis of different kinds of data, i.e., • an analysis of hydrological and hydromorphological data was collected during the field surveys; • development of conditional habitat suitability criteria according to the data of field surveys and fish monitoring; • evaluation of reference and altered (under influence of HPP) conditions according to historical observations and method of analogy; • simulation of the projections of pilot rivers using created HBV hydrological models that required the observed hydrometeorological data for calibration and validation, and output of ensemble of regional climate models according to different RCP climate scenarios for the discharge projections. The fourth level involved processing the mesoscale habitat modeling of the current situation using the MesoHABSIM model according to three main input components (river's hydromorphology, conditional habitat suitability criteria, and hydrological scenarios). In parallel, hydrological conditions under climate change were evaluated by combining the observed hydrological data and projected ones according to different RCP scenarios. Finally, the processing included the mesoscale habitat modeling in the future according to different hydrological scenarios. In the fifth step, the changes of fish habitat indices were estimated under both hydropower plants and climate change pressure.

Meso-Habitat Modelling
The physical habitat modeling system MesoHABSIM [45], was used to simulate spatial and temporal changes in aquatic habitat availability for different fish species under the influence of HPP at different climate change scenarios. MesoHABSIM aggregates data The fourth level involved processing the mesoscale habitat modeling of the current situation using the MesoHABSIM model according to three main input components (river's hydromorphology, conditional habitat suitability criteria, and hydrological scenarios). In parallel, hydrological conditions under climate change were evaluated by combining the observed hydrological data and projected ones according to different RCP scenarios. Finally, the processing included the mesoscale habitat modeling in the future according to different hydrological scenarios. In the fifth step, the changes of fish habitat indices were estimated under both hydropower plants and climate change pressure.

Meso-Habitat Modelling
The physical habitat modeling system MesoHABSIM [45], was used to simulate spatial and temporal changes in aquatic habitat availability for different fish species under the influence of HPP at different climate change scenarios. MesoHABSIM aggregates data of hydromorphological units (HMUs) and biological models, thus describing the characteristics and distribution of mesoscale habitats (i.e., glides, pools, riffles, rapids, and aquatic vegetation) [46], and linking them to the specific requirements of aquatic organisms in terms of a suitable range of depth, flow velocity, substrate composition, presence of covers and shelters, and other attributes that create suitable or optimal conditions for the species [47]. The habitat requirements for a species (biological model) are defined as habitat suitability criteria, developed using an inductive or deductive approach [48]. The habitat model can quantify the extension of the available habitat for target species by assigning a rank of suitability to each meso-habitat, according to the habitat suitability criteria for considered fish. For the MesoHABSIM application, time series of mean daily water discharge in natural and modified conditions for the year with all possible amplitude of target discharges were created to describe the characteristics of hydromorphological units in all possible hydrological conditions. Physical habitat parameters (depth, flow velocity, substrate composition, presence of shelters such as woody debris, boulders, aquatic vegetation, etc.) were collected during several field surveys at different flow conditions for each HMU identified within the selected river reaches. Habitat suitability criteria are developed to link the species distribution to these physical habitat parameters. Finally, to quantify the habitat evolution among space and time for target species, the SimStream plugin of QGIS was used [49].

Evaluation of Hydrological Characteristics
Commonly, the evaluation of river discharge during the low-flow period is being carried out using minimum values of the fixed averaging period. The hydrological characteristic of the average minimum discharge of 30 continuous days (Q 30 ) is used to characterize the low-flow conditions in Lithuania. Even environmental flow is estimated as 80% or 95% probability of Q 30 values according to national legislation. Hence, for the evaluation of hydromorphological and hydrological changes and later for fish habitat alterations, the four different low-flow situations (minimum of Q 30 , average of Q 30 , maximum of Q 30 , and annual mean) were chosen. These values were essential for the reasonable evaluation as all selected values covered all possible ranges of river discharge fluctuations during the warm period (May-October). First of all, a minimum of Q 30 and annual mean discharge (Q annual_mean ) were calculated for each year from the long-term dataset. According to the estimated values, the absolute values of Q 30_min , Q 30_ave , Q 30_max , and Q annual_mean were selected from the analyzed periods as target discharges for the measurements of field surveys and the range of possible discharge fluctuations during the low-flow period for the proper evaluation of fish habitats.
The hydrological scenarios for the reference (natural, without HPP) and altered (under influence of HPP) conditions were created for each river to estimate the impact of HPPs on fish habitats. For the Verknė River, the reference hydrograph was already available, since WGS is located upstream of the HPP. Therefore, the altered hydrograph was created according to the rules of HPP reservoir exploitation and energy production. Meanwhile, on the Širvinta, Šešupė, and Bartuva rivers, the water gauging stations (WGSs) are located downstream of HPPs; therefore, the altered hydrographs were already present because the observations of those WGSs were affected by HPPs operation. However, the reference hydrograph still needed to be prepared. For this task, the analogy method was chosen. This method establishes some logical correlation between the flow of two similar basins, one of which has not been disturbed or altered by the external intervention [50]. Two additional WGSs of Dubysa-Lyduvėnai (for the Širvinta and Šešupė rivers) and Minija-Kartena (for the Bartuva River) were selected for the restoration of natural conditions of river discharge at selected river stretch. The selection of a river analog was based on a similar catchment area, similarity in physico-geographical and hydrometeorological characteristics, and absence of anthropogenic structures, which could interrupt the continuity of the river, e.g., dams [39]. The regression equations between analyzed rivers and river-analog were created using daily discharge data before the date of HPPs construction ( Table 1). The mentioned equations and allocation of the annual water volume were used to reconstruct the reference hydrographs after the construction of HPPs. For each river, a representative year with the full range of possible low-flow fluctuations (from Q 30_min up to Q annual_mean ) during the warm season was selected to evaluate changes in indices of fish habitats.

Field Surveys of Hydromorphology
During the field works, the hydromorphology of the selected river reaches was investigated below the considered HPPs. To properly record the spatial heterogeneity of the HMU within each reach, we surveyed river stretches with a length of 10-20 times the river's widths. All the hydromorphological units (HMUs, also indicated as meso-habitats) identified within each river stretch were systematically mapped and digitalized in the shp format during the surveys in the season of 2020/2021 for the Verknė, Širvinta, and Šešupė rivers, and in 2017 for the Bartuva River. The measurements were done as close as it is possible to the warm period (May-October). However, due to the absence of higher discharges (Q 30_max and Q annual_mean ), some situations were measured in November and April. HMUs acquisition was achieved by means of a laser range-finder connected to a field computer via Bluetooth. For the spatial data collection and storage, Mapstream, a QGIS plugin designed for field data acquisition in the framework of the MesoHABSIM approach, was used. The definition of the HMUs was achieved following the classification proposed by [46], where glide, pool, riffle, rapid, and aquatic vegetation types are described. To assess the habitat evolution among different flows, river surveys were performed as close as possible to three different discharge conditions during the low flow period (Q 30_min , Q 30_ave , and Q 30_max ) and to one about the mean annual discharge. Physical habitat parameters such as depth and flow velocity within each HMU were collected by using the Valeport model 801 electromagnetic flow meter, mounted on a wading rod. In particular, in each HMU, 10 measurements of river depth and flow velocity were made, and, in the same points, substrate composition was assessed in terms of prevalent class. Finally, the presence/absence of covers and shelters for fish was assessed for each HMU.

Projection of River Runoff
The HBV hydrological model was applied to project river discharge under different RCP scenarios. The obtained alterations allow estimating the changes of fish habitat indices in the 21st century for the selected rivers (Verknė, Širvinta, Šešupė, and Bartuva) with the different hydrological conditions. The HBV model is a semi-distributed conceptual model that was developed at the Swedish Meteorological and Hydrological Institute (SMHI). This model is based on the water balance equation [51]: where P-precipitation, E-evaporation, Q-discharge, SM-soil moisture SP-snow pack, UZ-upper groundwater zone, LZ-lower groundwater zone, V-lake or dam volume. Model computations were performed in three steps. These steps include estimating the amount of precipitation that falls to the ground, estimating the slope runoff and the runoff in the watercourse as well as the runoff transformation. Calibration of the created models was performed using 16 basic parameters. The calibration procedure has to be processed until the correlation coefficient r becomes the highest and the relative difference (RE, %) between the measured and modeled discharge becomes the smallest. Calibration of the hydrological models was made for the period of 1986-1995 and validation for the period of 1996-2005. More about the application of the HBV model for the projections of Lithuanian river runoff is presented in many other studies [41,42,52,53]. The r for the calibration period varied from 0.70 to 0.83 and for validation periods-from 0.66 to 0.78 for the selected rivers. The RE in the calibration period varied from −7.3 to −10.4, and in the validation period varied from 7.7 to 13.3 offset each other ( Table 2). High r and small RE values allowed the use of hydrological models to project the runoff of selected rivers under climate change conditions. Despite the precise size of the RCM grid cells (11 × 11 km) of regional climate models, in some cases, their output gains systematic errors due to weak representation of distinctive local climatic conditions. Accordingly, the data of daily air temperature and precipitation amount was performed and adjusted to the Lithuanian conditions according to the quantile mapping method [54,55]: where MP Obs -observed meteorological parameter, MP RCM Ref -meteorological parameter simulated by the regional climate model for the reference period, ECDF Obs -empirical cumulative distribution function for the observed period, ECDF RCM Ref -empirical cumulative distribution function for the reference period of the regional climate model, MP RCM Fut -meteorological parameter simulated by the regional climate model for the future period. Historical observations from meteorological stations for the period 1986-2005 were used to adjust the RCM output of the particular grid cells.
To create future hydrographs under climate change conditions, the differences in days of target discharges (Q 30_min , Q 30_ave , Q 30_max , and Q annual_mean ), as well as the change of the target discharges, were calculated between modeled discharge of RCM future and RCM reference for different RCP scenarios. The estimated differences were incorporated in the altered hydrograph of the current situation.

Fish Habitat Models and Impact Assessment
Conditional habitat suitability criteria (CHSC) have previously been developed for 4 fish species: roach, dace, schneider, and vimba [39]. For the remaining species, within the current study, new CHSC were developed (Appendix A, Table A1) and validated using the same method [39]. Validation was carried out in 80 additionally sampled HMUs, with the exception of juvenile trout (72 HMUs) and juvenile salmon (22 HMUs), since not all river sections selected for validation for these species are accessible and/or habitable.
The suitable habitat area for the species (expressed in m 2 ) was modeled at reference and current altered conditions (under HPPs functioning) in a representative year and at altered conditions of the same year in the near and far future under different scenarios of climate change. The current impact of HPPs on habitat availability was assessed based on the comparison of the modeled available habitat area at reference conditions (baseline reference conditions) and under HPPs functioning. The combined impact of HPPs and climate change was assessed based on the comparison of the modeled available habitat area at the current altered conditions (baseline altered conditions) and altered conditions under different climate change scenarios. The flow value that exceeded 97% of the time at baseline conditions (Q 97 ) [56] and the corresponding area of species habitat (the minimum threshold area) were used as common denominators. Deviation of spatial habitat availability (hereafter, the index of spatial habitat availability; ISH) was quantified as the ratio between available habitat area in baseline and altered conditions. Deviation of temporal availability of suitable habitats was assessed based on the relative increase in the cumulative continuous duration of days when the area of the habitat falls below the minimum threshold values (hereafter, the stress days alteration; SDA), normalized between 0 and 1 by using the index of temporal habitat availability (hereafter, ITH). The concept and calculation of ISH, SDA, and ITH are described in more detail in [56][57][58]. The relative habitat area (% of river channel; S%) suitable for the modeled fish species and the relative habitat area at a discharge of Q 97 (SQ 97 %) at baseline reference conditions were also calculated. Only species with S% > 5 and SQ 97 % > 2.5 were used to assess the current and future impact of HPPs functioning in considered rivers.

Hydrological and Hydromorphological Changes
According to long-term datasets, the target discharges were calculated for each of the selected case studies sites (Table 3, column Hist). These discharges covered the potential range of discharge fluctuations in a particular river during the low-flow period of the warm season. The obtained results disclosed regional differences of the selected rivers. The calculated discharges of the Verknė and Širvinta rivers showed the smallest distribution between the values. These rivers have the strongest groundwater feeding compared to the other two. Therefore, this regularity causes more even discharges during the low-flow period. Meanwhile, the Šešupė and Bartuva rivers were distinguished by relatively low values of low-flow discharges compared to the annual mean discharge. The wide amplitude of target discharges confirmed the dependency of the mentioned rivers from surface runoff. The mapping of hydromorphological units (HMUs) was carried out as close as it is possible to obtain discharges (Table 3, column Actual). The differences in the measured discharge situations originated depending on the natural variability and anthropogenic regulation-HPP operation. The latter determined the absence of the Q 30_min situation in the Šešupė River. However, the anthropogenic activity allowed detecting extremely low values due to the lower boundary of the hydropeaking in the Širvinta and Šešupė rivers. Moreover, five different discharge situations were measured in the Širvinta River. The records of lower boundaries and additional data provided an opportunity to evaluate a wider range of low-flow discharges and their impact on fish habitats in more detail. The alterations of HMUs were evaluated depending on the four different discharges in the selected river stretch (Figure 3). For the Verknė River, most of the mapped HMUs were identified as glide HMU during low discharge conditions (Q 30_min and Q 30_ave ). Several units of pools and rapids were detected. At higher discharges (Q 30_max and Q annual_mean ), some glides transformed to riffle and rapid. A slightly different HMUs composition was obtained for the Širvinta River, where riffle units were found at low discharges. Later, they transformed into rapids due to an increase in river discharge. The differences of HMUs during the low-flow average and low-flow maximum situations were almost unchanged, only at annual mean discharge some pools move to glide. The most dramatic differences were observed in the Šešupė River because a huge loss of the HMU polygons was recorded at extremely low discharge. Some of the survived polygons were barely 20 cm in depth; therefore, aquatic vegetation became very dense and formed an HMU of aquatic vegetation. At higher discharges, most of the studied reaches were covered by glides and rapids at the river banks; however, during the annual mean discharge, the rapids were located in the middle of the river bed. In the Bartuva River, many small HMUs were found at discharges of low-flow minimum and low-flow maximum. Some units transformed from riffle to rapid or merged to rapids due to an increase in discharge. The HMUs of selected stretch became homogenous at the discharge of annual mean because all units were indicated as glides, only river depths and flow velocities defined the spatial differences and boundaries of the polygons.
transformed from riffle to rapid or merged to rapids due to an increase in discharge. The HMUs of selected stretch became homogenous at the discharge of annual mean because all units were indicated as glides, only river depths and flow velocities defined the spatial differences and boundaries of the polygons.
Despite the obtained changes in HMUs composition with the respect to the different discharge conditions, the alterations in the hydraulic features (river depth and flow velocity) were the greatest. These parameters were one of the most important characteristics that caused the changes in indices of fish habitats due to discharge fluctuations.  Despite the obtained changes in HMUs composition with the respect to the different discharge conditions, the alterations in the hydraulic features (river depth and flow velocity) were the greatest. These parameters were one of the most important characteristics that caused the changes in indices of fish habitats due to discharge fluctuations.

Projections of Rivers Discharge and Hydrological Scenarios
The projections of discharge of the Verknė, Širvinta, Šešupė, and Bartuva rivers were made for the near (2021-2040) and far future (2081-2100) according to three RCP scenarios (2.6, 4.5, and 8.5). The most attention was paid to target discharges of the warm season low-flow (Q 30_min -Q annual_mean ) because the variability of them determined hydromorphological and hydraulic conditions in defined HMUs, while the fish habitat alterations were closely related to hydromorphological conditions. The primary analysis consisted of the future changes in the average number of days with target discharges under climate change compared to the reference period (Table 4). According to projected values, the increase in days of Q 30_min and Q 30_ave was obtained. In the near future, the highest increase was found according to RCP2.6 and RCP8.5 in the Šešupė River. Here, the average number of days of Q 30_min increased by 8 and 9 days, respectively in comparison to the reference period. Similar growth in days Q 30_min was estimated in Verknė River according to RCP 8.5. In the Širvinta and Bartuva rivers, days with the lowest discharges were almost unchanged because the number of days fluctuated between the decrease and the increase in three days. The situation changed in the far future when clear an increase in days of Q 30_min was estimated for all rivers. Especially, the differences between RCP scenarios gained consistent patterns. The smallest increase by 5-8 days was projected according to the RCP2.6 scenario, while the highest differences (15-32 days) compared to the reference period were found according to RCP8.5. The Šešupė River was obtained as the most affected river under climate change conditions because the projections showed the largest increase in days with Q 30_min for all RCP scenarios against the other rivers. An increase in days of Q 30_ave was found for all rivers except the Širvinta. The estimated values of Q 30_ave highlighted the similar distribution of alterations between RCP scenarios as in the case of Q 30_min . The obtained increase in the average number of days with the lowest discharges was caused by the overall reduction of rivers discharge. Consequently, the decrease in the average number of days of Q 30_max and Q annual_mean was projected. In terms of the changes in the quantity of river discharge, the alterations of projected discharges were also evaluated. The deviations for target discharges were evaluated for the same future periods and according to the same RCP scenarios as in the case of the number of days (Figure 4). The hugest changes were obtained for the Q 30_min and Q 30_ave situations compared to the reference period. The estimated deviations fluctuated between −60% and 40% for the Q 30_min , and in the range of ±20% for Q 30_ave . The discharges of annual mean and all other remaining (higher) discharges were the least affected. The selected rivers were mostly distinguished from each other by Q 30_min because in the near future the Q 30_min of Verknė River was almost unchanged. Whereas in the Šešupė River the projected changes showed up to −30% deviation for all RCP scenarios. For the far future, the decrease of Q 30_min was projected for Verknė, Širvinta, and Šešupė rivers. The climate change effect on mentioned rivers distributed evenly, i.e., the smallest decrease of Q 30_min was found in the Verknė river, moderate-in the Širvinta, and the highest decrease-in the Šešupė River according to RCP4.5 and RCP8.5. The Bartuva was unique compared to the other rives because, in the near future, the increase of Q 30_min and Q 30_ave , and the decrease of Q 30_ave were estimated. Quite similar deviations were projected for the far future; only Q 30_min slightly decreased (up to 10%) under the RCP4.5 and RCP8.5 scenarios. The decrease of Q 30_ave values was projected for the Bartuva River in the near and far future. The analyzed characteristics of change in days with target discharges and deviation in values of those discharges were important for the evaluation of fish habitats as ISH and ITH indices were related to the quantity and temporal availability of the target discharges. Hence, the changes in those discharges directly affected alterations of the mentioned indices.
crease of Q30_min was projected for Verknė, Širvinta, and Šešupė rivers. The climate change effect on mentioned rivers distributed evenly, i.e., the smallest decrease of Q30_min was found in the Verknė river, moderate-in the Širvinta, and the highest decrease-in the Šešupė River according to RCP4.5 and RCP8.5. The Bartuva was unique compared to the other rives because, in the near future, the increase of Q30_min and Q30_ave, and the decrease of Q30_ave were estimated. Quite similar deviations were projected for the far future; only Q30_min slightly decreased (up to 10%) under the RCP4.5 and RCP8.5 scenarios. The decrease of Q30_ave values was projected for the Bartuva River in the near and far future. The analyzed characteristics of change in days with target discharges and deviation in values of those discharges were important for the evaluation of fish habitats as ISH and ITH indices were related to the quantity and temporal availability of the target discharges. Hence, the changes in those discharges directly affected alterations of the mentioned indices.  Fish habitat evaluation with meso-habitat modeling required the datasets of daily discharge for the reference (natural) and altered conditions. Therefore, the different hydrological scenarios of daily discharges were systematized and prepared not only for the reference and altered conditions in the historical period but also under the selected climate change scenarios ( Figure 5). The year representing all amplitude of target discharges was selected for each case study river. The near and far future hydrographs were created according to previously obtained alterations of the projections of target discharges. Altered hydrographs of historical observations were used as a background and reference point evaluating climate change impact. The upper boundary of the results was displayed as the discharge of the annual mean because the MesoHABSIM model does not evaluate the values above the upper threshold. The results of hydrological scenarios clearly indicated a crucial impact of HPPs on the natural river flow. The hydropeaking was expressed instead of the natural variability in the Verknė, Šešupė and especially Širvinta rivers in the historical period. The decrease in the lower boundary of discharges was the most noticeable effect of HPP regulation compared to the reference conditions. For the Bartuva River, the impact of HPP on the river flow was the weakest, as the river itself had a natural variability due to rainfall events. The main effect was reflected on the low-flow of the warm season when the altered river discharge gained lower values under the redistributed runoff due to HPP operation. Despite the effects of the anthropogenic activities, the joint HPP and climate change influence was studied. The results indicated the Šešupė River as the most affected river because in the near future the ensemble of RCP scenarios projected a dramatic decrease of its discharge during the warm season. In the far future, the same ensemble projected even greater changes of negative signs. Meanwhile, the Verknė and Širvinta rivers were less affected; only in the far future, the combined effect emerged as a decrease

Fish Habitat Changes under the Influence of HPP and Climate Change
The modeling results show that, at present, the functioning of HPPs does not significantly affect the spatial availability of suitable habitats for fish species in all studied rivers (ISH values are in the range 0.87-1.0) ( Table 5). The impact on the temporal availability of habitat is greater, but the cumulative duration of events where the area of suitable habitat falls below the threshold area available for fish under natural conditions at Q97 (stress days alteration SDA) did not exceed 100% for most modeled species, with the exception of juvenile salmon and juvenile trout in the Verknė, schneider, minnow, and chub in the Šešupė, Širvinta, and Bartuva and gudgeon in the Bartuva (SDA up to 156%, ITH ≥ 0.55) ( Table 5).

Fish Habitat Changes under the Influence of HPP and Climate Change
The modeling results show that, at present, the functioning of HPPs does not significantly affect the spatial availability of suitable habitats for fish species in all studied rivers (ISH values are in the range 0.87-1.0) ( Table 5). The impact on the temporal availability of habitat is greater, but the cumulative duration of events where the area of suitable habitat falls below the threshold area available for fish under natural conditions at Q 97 (stress days alteration SDA) did not exceed 100% for most modeled species, with the exception of juvenile salmon and juvenile trout in the Verknė, schneider, minnow, and chub in the Šešupė, Širvinta, and Bartuva and gudgeon in the Bartuva (SDA up to 156%, ITH ≥ 0.55) ( Table 5).
Modeling future scenarios shows different patterns in the deviation of the spatial and temporal availability of habitats suitable for certain fish species in different rivers. Irrespective of the climate change scenario, the reduction in the spatial availability of the habitats due to the operation of HPP in the future, compared to the current altered conditions, will remain relatively insignificant (ISH > 0.90) for all fish species in the Verknė, Širvinta, and Bartuva rivers. A larger spatial reduction in suitable habitats is predicted only in the Šešupė River in the far future under the scenarios RCP4.5 and RCP8.5, where ISH values for all modeled fish species are in the range 0.89-0.78 (Table 6). But significant changes are foreseen in the temporal availability of habitats. In the Verknė River, according to climate change scenarios RCP4.5 and RCP8.5, a more than 2-fold increase in the duration of events (ITH ≤ 0.48) is predicted when the area of habitats suitable for 7 out of 10 modeled species temporarily drops below the current SQ 97 % in the far future. The number of stress days for chub, bullhead, vimba, and juvenile salmon is also expected to increase more than 1.5-fold (ITH < 0.57) in the near future at RCP8.0 and in the far future at RCP8.5. In the Širvinta River, the temporary change in the area of suitable habitat will remain negligible for all species at the scenarios RCP4.5 and RCP8.5 in the near future, but more than a half of the modeled fish species will be affected by the temporary reduction in the suitable habitat at the scenario RCP2.6 in the near future and RCP4.5 and RCP8.5 in the far future. At the scenario RCP2.6, a more than 3.5-fold increase in the number of stress days (ITH < 0.27) for dace and juvenile salmon is also predicted in the far future. In the Šešupė, 4 out of 6 species whose suitable habitat has been modeled in this river will be affected by a significant increase in the number of stress days at all climate change scenarios in both the near and far future periods. The most significant increase in the number of stress days will be in the far future at the scenarios RCP4.5 and RCP8.5. According to the last two climate change scenarios, the number of days, during which the area of suitable habitat for schneider, dace, minnow, and chub falls below the threshold area available under current conditions at Q 97 %, will increase by more than 5 times (ITH < 0.15). The increase in the number of stress days will remain insignificant only in the Bartuva River for all fish species, regardless of the climate change scenario. Table 5. The relative mean habitat area (% of river channel; S%) of the modeled fish species, the relative habitat area at a discharge of Q 97 (SQ 97 %) at natural conditions, and the average amount of habitat loss (index of spatial habitat availability; ISH), and the relative increase (in %) in the cumulative duration of stress days (stress days alteration SDA and corresponding index of temporal habitat availability ITH) when the area of habitat falls below SQ 97 % when HPPs function. Species with a modelled S% < 5 or SQ 97 % < 2.5 are excluded. Compared to the average ITH for all modeled fish species in a given river under a specific RCP scenario, the increase in the number of stress days when the area of habitat suitable for gudgeon and stone loach falls below the threshold compared to the baseline in the Širvinta and Šešupė rivers is relatively less than for other fish species under all climate change scenarios ( Figure 6). The same is true for the Venta River, with the exception of the stone loach at the scenario RCP8.5 in the far future. The impact of the increased stress days for minnow and juvenile trout in the Verknė and that of bullhead in the Širvinta will also be relatively smaller compared to other fish species. The most significant negative impact of HPP operation in the context of climate change will be on the temporary availability of habitats suitable for schneider, dace, chub, juvenile salmon, and vimba in all of the above-mentioned rivers, where these species should live in natural conditions. Only in the Bartuva River, the relative increase in the number of stress days is similar for all fish species, except for the schneider, whose ITH values deviate to a greater extent depending on the climate change scenario.

Discussion
The overall tendencies disclosed the high vulnerability of the Šešupė River from the Central hydrological region, where rivers are closely dependent on the surface feeding, especially on snowmelt. Climate change is likely to increase the number of winters without snow cover due to raising the temperature; therefore, the natural response to such change was felt even on the discharges of the warm period, especially on the lowest discharges. Meanwhile, the rivers (Verknė and Širvinta) from the South-eastern hydrological   ISH ITH ISH ITH ISH ITH ISH ITH ISH ITH ISH ITH ISH ITH ISH ITH ISH ITH ISH ITH ISH ITH ISH

Discussion
The overall tendencies disclosed the high vulnerability of the Šešupė River from the Central hydrological region, where rivers are closely dependent on the surface feeding, especially on snowmelt. Climate change is likely to increase the number of winters without snow cover due to raising the temperature; therefore, the natural response to such change was felt even on the discharges of the warm period, especially on the lowest discharges. Meanwhile, the rivers (Verknė and Širvinta) from the South-eastern hydrological region suffered less from the impact of low flow due to compensating effect of groundwater feeding. The projected amount of liquid precipitation under climate changes in some cases could increase despite the negative effect of temperature on snow events [42,52]. Accordingly, the Bartuva River from the Western hydrological region (rainfall-depended river) was almost unchanged, or the projected conditions revealed an average increase in the lowest discharges according to certain scenarios. The obtained regularities once more confirmed the importance of the natural origin and behavior of each river. This indicates that regional differences should be considered in the evaluation of projections under climate change.
The modeled spatial and particularly temporal availability of habitats suitable for different fish species also disclosed different patterns in the rivers from different hydrological regions. However, the comparison of the deviation of ITH values from the average ITH value for all fish species in a certain river suggests that small bottom dwellers, gudgeon, and stone loach are likely to become increasingly dominant in all studied rivers, except the Bartuva River. The relative increase in the cumulative duration of stress days when the area of habitat suitable for these fish species falls below SQ 97 % is less than for other modeled species that are common in Lithuanian rivers. The stress days alteration for minnow and juvenile trout in the Verknė, and that for bullhead in the Širvinta is also less than average for all fish species. Such patterns correspond to findings [59], who showed that the operation of small HPPs in Czech rivers caused the replacement of large-bodied fish species (adult brown trout, chub, dace, and grayling) towards small-bodied fish (juvenile trout, minnow, bullhead, stone loach, and gudgeon). Ovidio et al. [60], who studied the impact of a newly built hydropower plant and the corresponding increase in the duration of low flow, also reported a decrease in the biomass of European grayling (Thymallus thymallus) and brown trout and a shift in their age structure towards younger fish, while the number of stone loach and bullhead, on the contrary, increased significantly after several years of HPP operation.
Both the results of modeling of habitat availability under climate change conditions in rivers affected by hydropower plants and publications on empirical studies show that operation of HPP can be beneficial for small fish such as gudgeon, bullhead, or stone loach in competition with other species. Interspecies competition may also interact with climate change [61,62], while the impact of HPP leads to additional selection through different effects on the temporal availability of habitat for different fish species. Habitat heterogeneity is also a very important factor in determining the resilience of fish to changes in flow [63,64], therefore the impact of HPP on fish under climate change may be most pronounced in the rivers with lower habitat diversity.
In the future, along with an increase in the number of days with the minimum flow, an increase in water temperature is expected. This could overshadow the advantage of small-bodied stenothermal fish due to the relatively lower temporary loss of habitat during the operation of the hydropower plant. Trout, bullhead, and minnow, as well as schneider and salmon, are stenothermic fish [65,66]. An increase in temperature may outweigh the relative gains of small stenothermic fish, and further, enhance the negative effect of HPP on larger ones. Conversely, the rising temperature can mitigate HPP-induced temporary reduction of suitable habitat for bleak, roach, chub, and other cyprinid fish species that prefer or tolerate higher water temperatures. A combination of lower discharge with higher temperature has been shown to have a positive effect on the recruitment of some of these cyprinids [67].
The simulation results show that depending on the climate change scenario and the river region, the impact of the HPP on temporal habitat availability for fish will become even more significant. Already in the near future, this impact should be most noticeable in the HPP-affected Lithuanian rivers of the Central region, which discharge is dominated by surface runoff. In the groundwater-fed rivers of the South-eastern region, a significant increase in the impact of HPP on temporal availability of suitable habitats for certain fish species is more likely in the far future, while in the near future, the possibility of a significant impact is higher only at the RCP2.6 climate change scenario. The exception is the rivers in the Western region (represented by the Bartuva River in this study), where small changes in the impact of HPP on fish habitat are predicted. Only the temporary availability of habitat suitable for a schneider may deviate depending on the climate change scenario.  Institutional Review Board Statement: Ethical review and approval were waived for this study. Fish sampling was performed according to the European Standard CEN EN 14011:2003 "Water quality -Sampling of fish with electricity". All fish caught remained healthy and alive and were released immediately after species identification and counting. The permit to sample fish with electric fishing gear was issued by the Lithuanian Environmental Protection Agency (permit No. 033, issued on July 18, 2017 m. and permit No. 026, issued on June 16, 2020).

Data Availability Statement:
The data that support the findings of this study are available from the paper authors upon reasonable request.