System Dynamics Modeling of Water Level Variations of Lake Issyk-Kul, Kyrgyzstan

Abstract

Lake Issyk-Kul is an important endorheic lake in arid Central Asia. Climate change, anthropogenic water consumption and a complex basin hydrology with interlocked driving forces have led to a high variability of the water balance and an overall trend of decreasing lake water levels. The main objective of this study was to investigate these main driving forces and their interactions with the lake’s water level. Hydro-meteorological and socioeconomic data from 1980 to 2012 were used for a dynamic simulation model, based on the system dynamics (SD) method. After the model calibration and validation with historical data, the model provides accurate simulation results of the water level of Lake Issyk-Kul. The main factors impacting the lake’s water level were evaluated via sensitivity analysis and water resource scenarios. Results based on the sensitivity analysis indicated that socio-hydrologic factors had different influences on the lake water level change, with the main influence coming from the water inflow dynamic, namely, the increasing and decreasing water withdrawal from lake tributaries. Land use changes, population increase, and water demand decrease were also important factors for the lake water level variations. Results of four scenario analyses demonstrated that changes in the water cycle components as evaporation and precipitation and the variability of river runoff into the lake are essential parameters for the dynamic of the lake water level. In the future, this SD model can help to better manage basins with water availability uncertainties and can guide policymakers to take necessary measures to restore lake basin ecosystems.

1. Introduction

Water is considered the most precious natural resource [1], and its availability is directly threatened by human activities [2]. In arid regions, water scarcity is already today of the highest significance while ongoing socioeconomic growth and the climate change will only intensify the challenges related to managing the increasingly limited water resources. Climate change leads to an increase of the air temperature and more variable rainfall regimes, with severe consequences for the frequency and magnitude of droughts and flood events, and an accelerated meltdown of glaciers, which can increase the river runoff in the short term but ultimately alters the discharge regimes in the long term, reducing the amount of available water [3,4]. Combined with projected population growth, expanding irrigation schemes, increased evapotranspiration rates, and the growing water demand for domestic and industrial purposes, these climatic changes will have a significant impact on the water balance in the near future [5]. This thematic issue is particularly serious in Central Asia [6], a region characterized by intense water-related issues and heavily affected by the climate change [7,8,9].

Most of Central Asia is determined by an arid and semiarid climate. Therefore, areas close to shallow groundwater, rivers, and lakes are characterized by unique water-dependent ecosystems and human societies, which are closely related to the limited water resources and have been developed for thousands of years [10]. At the beginning of the 21st century, global climate change, population growth, high agricultural water use, rising levels of pollution, river flow regulation and dropping lake levels exert multiple pressures on the Central Asian water resources [11,12,13]. It is evident that in many lake catchments the agriculture development is associated with increasing water withdrawal from lake tributaries, which causes serious lake water level drops [14,15,16].

Endorheic lakes located in arid and semi-arid environments are highly sensitive to climate oscillation since minor changes in their environmental conditions can cause dramatic changes of the lake level [17]. The sensitivity of these lake is even higher when they are located in high-altitude mountain environments [18]. In addition, many researchers indicated that the notable climate warming and decrease in precipitation since the 19th century is associated with an abrupt lake water level drop in arid regions worldwide [19,20,21]. Lake Issyk-Kul is a prime example of these endorheic mountain lakes in arid Central Asia, and the variations of its water level reflect the global climate processes very well [22].

Lake Issyk-Kul is a closed lake, located in the Tian-Shan mountain belt within the Republic of Kyrgyzstan, in arid Central Asia [23]. From 1973 to 1999, the average monthly temperature in the lake basin increased, resulting in an acceleration of glaciers melting and an increased river runoff during the summer period. However, the lake level has decreased by about 13 m since in the mid-19th century, while 3 m of lake level decrease has been recorded from 1927 to 1966 [21]. In recent years, Lake Issyk-Kul became a subject of special concern because of the continuous water level decrease and its contamination related to human activities [24]. The variations of the lake water level have impacts on the agriculture, the local population, economic activities and the environment [25,26,27]. Studies conducted on the Issyk-Kul during the last decades covered topics ranging from geological characteristics and earthquakes in the Northern Tian-Shan mountain range [28]; the reconstruction of the paleoclimate based on sediment cores and isotopes [29,30]; the lake basin glaciation and the impacts of the climate change [31]; the physicochemical characterization of the basins water resources in general [32,33] and the influence of uranium mining in particular [34,35], to the potential for a biodiversity-friendly fishing industry in the lake [36]. Specific components of the Issyk-Kul water balance have also been studied by several projects. The main focus of those studies has been on the groundwater resources and the deep-water renewal as a part of the Issyk-Kul water balance [37,38,39,40,41], while only very few analyses of the lake level variability have been published [16,21,42,43,44]. These previous studies dealing with the water balance are too limited in their focus, as they only analyze a single aspect of it (groundwater, climate change impact on glaciers, or the physicochemical properties of the lake’s water) of the complex hydrological lake catchment interactions. They neglect equally important aspects like the land cover and its dynamic, human activities or the variability in precipitation discharge. As well, that the few studies focusing on the lake level variations, as an integral of the water balance, cover a temporal scale that is more related to paleoclimatic research and the influence of the last glaciation then on the present day water resources problems.

In order to close this knowledge gap and to allow a better understanding of the lake level variations and its main driving forces a complex, holistic and dynamic simulation model of the lake and its whole catchment is needed. However, these requirements seem difficult to be met as the system is highly complex and its individual parameters are controlled by numerous cause and effect relations, with feedback mechanisms and nonlinear characteristics [45,46,47]. One method suitable for this ambitious task is the system dynamics (SD) model, which was originally developed by Forrester in 1961 [19], and is an approach for understanding the interactions among driving factors and interconnected sub-systems that drive the dynamic behavior of a system [48,49]. Over the years, a number of SD models have been developed for water balance simulation and have been used to evaluate various water-related solutions [50,51,52], such as water resource planning models [53,54,55,56], hydrologic extremes models [57], agriculture water management models [58,59], and water balance models, which have been developed to test water-related and environmental issues in developing countries where the data availability is lacking [60]. With this background, the SD model satisfies the requirements for a complex analysis of the Issyk-Kul water level fluctuations and its driving factors.

The novelty of the present study is the holistic approach which incorporates the whole lake catchment and a wide range of parameters influencing the water balance, including land use, the population, the economy, water supply and demand, and how these parameters change over time and interact with each other. To achieve this goal, a dynamic water balance model for the period 1980 to 2012 was developed in a first step, and the relations between the relevant factors were formulated based on mathematical equations. The second step was the calibration and validation of the model using historical data. The third and final step was to evaluate the main factors impacting the lake water level via a sensitivity analysis and various water resource scenarios.

2. Model Structure

2.1. System Dynamics Model

The event-oriented view of the world and linear thinking cannot sufficiently address complex issues. The SD model, based on the nonlinear causal thinking, has been developed to expand the dynamic simulation model [60]. In this study, based on the SD principles and a series of casual loops and mathematical equations, a novel stock-flow diagram has been developed for the Issyk-Kul basin, considering the hydrological and socioeconomic drivers of the lake water level change and their interactions (Figure 1). This model uses the concepts of water balance and rainfall-runoff transformation and estimates the water availability for the uses of agricultural, domestic, and industrial demands. Based on land use changes, population growth, and industrial productivity progression, the model generates water demand data for agricultural, domestic, and industrial use, respectively. The SD model also improves the understanding of the internal working of a complex socio-hydrological system with its various direct and indirect interconnections.

The SD model is divided into seven major sectors: water volume of the lake; groundwater volume; land use changes; population, and; water demands for agricultural, domestic and industrial use. Each sector defines an important component of the system. The relationships and the feedback of the stock-flow diagram are developed in the Vensim DSS^® software [61], as shown in Figure 1 to characterize the system processes. The “+” marks close to the arrows show an increase of the variable at the head of the arrow, while the “−” marks indicate a decrease of that variable. The time horizon for the simulation model extends from 1980 to 2040, and the model operates on an annual time step. The period from 1980 to 2012 is used for the model calibration and the remaining period (2013–2040) for the simulation of possible future developments.

2.2. Model Equations

In the equations of the SD model, the variables are either stocks or flows. Stocks are accumulations and are used to characterize the state of the system and create the information upon which decisions are made. Flow variables are used to define rates which can change the stock variables [62]. Stocks accumulate the difference between its inflow and outflow. Therefore, a stock with one inflow and one outflow is formulated as follows:

$$Stock (t) = {{\displaystyle \int }}_{t_0}^t \left[Inflow \left(t\right) - Outflow \left(t\right)\right]\mathrm{d}t + Stock (t_0)$$

(1)

where Stock(t) = Amount of stock at time t, Inflow(t) = Inflow at time t, Outflow(t) = Outflow at time t, and Stock(t₀) = Stock in time t₀. According to the Equation (1), other mathematical equations of the major processes of inflow and outflow represented in the following sections from Section 2.2.1, Section 2.2.2, Section 2.2.3, Section 2.2.4, Section 2.2.5.

2.2.1. Lake Water Volume

The Issyk-Kul water volume (storage) is the main stock of the model. Annual changes of the lake’s water volume are controlled by the in- and outflow of water. The inflow to the lake includes precipitation, surface and groundwater inflow. As the Issyk-Kul is an endorheic lake the outflow is limited to evaporation. The amount of water in the Issyk-Kul at a time t is calculated as follows:

$$V \left(t\right) = V (t_0) + {{\displaystyle \int }}_{t_0}^t \left[ GW \left(t\right) + (SW (t) - W (t\left)\right) + P \left(t\right) - E \left(t\right)\right]\mathrm{d}t,$$

(2)

where V(t) = Water volume of the lake at time t, V(t₀) = Water volume of the lake at the time t₀, GW(t) = Groundwater into the lake at time t, SW(t) = Surface water flow at time t, P(t) = Precipitation at time t, and E(t) = Evaporation from the lake surface at time t. W(t) = Water consumption from lake tributaries at time t, which is divided into Q₁, Q₂, and Q₃—is the amount of water used for domestic, agricultural and industrial, respectively. Surface water inflow = ((SW) − W(t)).

2.2.2. Groundwater

Groundwater inflow is an important component in the lake water budget and thus was first considered as a separate component of the Issyk-Kul water balance by Kaplinsky and Timchenko [63], who showed that the groundwater component provides about 30% of the income. The groundwater aquifers are located in permeable coarse Meso-Cenozoic sediments ((cemented) pebbles, sand, and clay as well as siltstone conglomerates) of up to 5 km in depth. The youngest groundwater layer can be found in up to 300 m-thick Quaternary sediments (consolidated coarse clastics) with a high permeability (0.04–4.0 m/h) and a total storage capacity of 58 km³. This Quaternary groundwater discharges into the lake through littoral and submarine vents and provides most of the groundwater inflow (>70%), though deeper groundwater aquifers (e.g., the Pliocene-Quaternary Sharpyldak or the Paleogene-Neogene Issyk-Kul aquifer) also contributes to the lake’s water balance [41]. The special topography and geological structure of the lake Issyk-Kul basin, as well as the groundwater gradient maps, show no evidence for groundwater outflow from the lake and its adjacent aquifers in the region [41]. Therefore, only the groundwater inflow of into the lake was considered in this study. The groundwater inflow of the lake was calculated with this equation:

$$GW \left(t\right) = GW(t_0) + {{\displaystyle \int }}_{t_0}^t \left[{R \left(t\right) - W \left(t\right) - GW}_L)\right]\mathrm{d}t$$

(3)

where GW(t) = Groundwater into the lake at time t, GW(t₀) = Groundwater into the lake at time t₀, R(t) = Recharge at time t, W(t) = Water demand/consumption at time t (it is the amount of water used for domestic, agricultural and industrial), and GW_L = Groundwater losses at time t (evaporation from the groundwater surface [41]).

The recharge primarily comes from the infiltration from surface water bodies, recharge from domestic, agricultural, and industrial water use, and percolation of precipitation. The recharge into the groundwater is calculated as:

$$R \left(t\right) = R (t_0) + {{\displaystyle \int }}_{t_0}^t \left[C_I \times P \left(t\right) + ( R_1 + R_2 + R_3) + I_{\mathrm{swb}}\right]\mathrm{d}t$$

(4)

where R(t) = Recharge at time t, R(t₀) = Recharge at time t₀, C_I = Coefficient of infiltration, P (t) = Precipitation at time t, R₁ = Recharge from domestic water demand, R₂ = Recharge from agricultural water demand, R₃ = Recharge from industrial water demand, I_swb = Infiltration from surface water bodies.

2.2.3. Water Demand

The GIP is a measure of the productivity of the population and hence is associated with the water consumption. In this model, water demand includes the domestic water demand, the agricultural water demand, and the industrial water demand. These can be calculated as:

$$T_{\mathrm{wd}}{ = D}_{\mathrm{wd}}{ + A}_{\mathrm{wd}} + I_{\mathrm{wd}}$$

(5)

where T_wd = Total water demand (in million m³/year), D_wd = Domestic water demand (in million m³/year), A_wd = Agricultural water demand (in million m³/year), and I_wd = Industrial water demand (in million m³/year).

The domestic water demand is the quantity of water used in people’s daily life. It can be estimated by the size of the population and the water demand per capita per day. The domestic water demand was formulated as follows:

$$D_{\mathrm{wd}}{ = p\left(t\right) \times D}_{\mathrm{pwd}} \times 365/{10}^3/{10}^8$$

(6)

where p (t) = Population, and D_pwd = Population water demand per capita (L/person·day).

The agricultural water demand is the water used in irrigation. It can be calculated by the agricultural area and the irrigation water demand per hectare calculated as follows:

$$A_{\mathrm{wd}}{ = A}_{\mathrm{pwt}}{ \times L}_{\mathrm{agr}}/{10}^8$$

(7)

where A_pwd = Irrigation water demand per hectare in cubic meters (m³/ha), and L_agr = Agricultural land area (ha).

The industrial water demand includes the water used in the process of production. Water use in industry varies significantly, the gross industrial product was used to calculate the industrial water demand. It is calculated as the gross industrial product multiplied by the average water demand for one unit of industrial product, formulated as:

$$I_{\mathrm{wd}}{ = I}_{\mathrm{pwd}} \times GIP/{10}^8$$

(8)

where I_pwd = Water demand for the industrial production (Kyrgyz Som (Ks)/m³), and GIP = Gross Industrial product (Ks).

2.2.4. Population

The inflows to the population stock consist of births and immigration; the outflows consist of deaths and emigration. The population and net population change in the lake basin are formulated as:

$${p \left(t\right) = p (t}_0{) + p (t}_0) \times (B\mathrm{r} - D\mathrm{r} + I\mathrm{r} - E\mathrm{r})$$

(9)

$$np \left(t\right) = p_{\mathrm{g}} \times p \left(t\right)$$

(10)

where p(t) = Population at time t, p(t₀) = Population at time t₀, Br = Birth rate, Dr = death rate, Ir = Immigration rate, Er = Emigration rate, np(t) = net population change, and p_g = population growth rate.

2.2.5. Land Use Changes

Land use patterns in the study area are divided into three categories-agricultural, natural, and urban land. Land use changes are a complex process, managed by several factors including population growth, economics, and legislation. In the SD model, it is not possible to foresee and simulate the policy changes, but the population growth and the resulting land use dynamics can be simulated. For this, it was assumed that, based on the population dynamics, natural land is converted into agricultural land, and natural land and agricultural land are both transformed into urban land. The urban land was determined by cumulating the changes in the agricultural land and natural land over time, calculated as follows:

$$L_{\mathrm{urb}}{\left(t\right) = L}_{\mathrm{urb}}(t_0) + {{\displaystyle \int }}_{t_0}^t \left[R_{\mathrm{agr}}{ + R}_{\mathrm{nat}}\right]\mathrm{d}t$$

(11)

$$L_{\mathrm{agr}}{\left(t\right) = L}_{\mathrm{agr}}{\left(t_0\right) + C}_{\mathrm{agr}} - {{\displaystyle \int }}_{t_0}^t{R}_{\mathrm{agr}}\mathrm{d}t$$

(12)

$$L_{\mathrm{nat}}{\left(t\right) = L}_{\mathrm{nat}}(t_0) - {{\displaystyle \int }}_{t_0}^t{R}_{\mathrm{nat}}\mathrm{d}t$$

(13)

The urbanization rates for agricultural land and natural land were used to estimate the total land use change into urban land. The urbanization rate is a function of land patterns, which is computed by the population growth and urbanization factors of natural land and agricultural land. These parameters were determined by analyzing the change in the total population and land use in the lake basin. The urbanization rates for natural and agricultural land in the lake basin were calculated using:

$$R_{\mathrm{nat}}{\left(t\right) = F}_{\mathrm{nat}} \times np\left(t\right)$$

(14)

$$R_{\mathrm{agr}}{\left(t\right) = F}_{\mathrm{agr}} \times np(t)$$

(15)

where L_agr = Agricultural land (ha), L_urb = Urban land (ha), L_nat = Natural land (ha), C_agr = Natural land converted into Agricultural land (ha), R_agr = Urbanization rates for agricultural land (ha/year), R_nat = Urbanization rates for natural land (ha/year), F_agr = Urbanization factor of agricultural land (ha/person), and F_nat = Urbanization factor of natural land (ha/person).

3. Materials and Methodology

3.1. Study Area

The Issyk-Kul (Figure 2) is an endorheic mountain lake, located on the northern slopes of the Tian-Shan mountains, Republic of Kyrgyzstan, in arid Central Asia (42°25′ N, 77°15′ E). It is situated at an altitude 1607 m above sea level, covering an area of 6236 km². The maximum depth of the lake is 668 m and the total water volume is 1736 km³. It is surrounded by high mountain ranges: the Kungey Ala-Too Range in the north with the highest peaks reaching 4770 m, and the Teskey Ala-Too Range in the south with peaks exceeding 5200 m [16]. More than 118 rivers feed into the lake. The main tributaries are the Pzhergalan River and the Tyup River, and their runoff regime is dominated by meltwater from glaciers and snow [23,64]. At present, the lake has no outlet. Up to the late Pleistocene, the Chu River was flowing into and out of the lake at its western extremity. As a result of a tectonic shift, the river changed its course less than five kilometers before reaching the lake and now flows in southwestern direction into the Boom Canyon and into the Naryn region [65]. Hence, the water budget of the lake is controlled by inflow into the lake and includes precipitation, surface water inflow, and groundwater inflow, while the outflow from the lake is limited to the evaporation.

3.2. Model Inputs

This study used annual precipitation data (1980–2012) collected from four meteorological stations (Balykchy, Karakol, Kyzyl-Suu, and Cholpan-Ata) located around the lake. The average of these stations was considered as the precipitation data inflow into the lake. Annual surface runoff from 16 hydrometric stations (Ak-Terek, Ak-sai, Ton, Tossor, Tamga, Dzhuuku, Chong-Kyzyl-Su, DZhety-Oguz, Karakol, Pzhergalan, Tyup, Chong-Uryukty, Ak-Suu, Chong-Ak-Suu, Cholpon-Ata, and Chong-Koi-Suu, see Figure 2), which are the main rivers leading up to the Issyk-Kul were used to determine the surface water flow. For the evaporation term, estimates of the mean annual energy budget at the water surface and evaporation from the lake surface were obtained from Salamat et al. [16]. The lake volume-area and volume-depth curve data were extracted from Shabunin and Shabunin (Figure 3) [22]. The groundwater inflow contributes about 30% of the water balance income [63], but no research has been done on how much groundwater actually flows into the lake. Therefore, in this paper, according to the Equations (3) and (4), the groundwater inflow is calculated. The resulting total groundwater inflow was then fitted, using STELLA’s graphical functions to estimate the unknown relevant factors of the water balance, like the coefficient of infiltration and water losses from groundwater. There is a small difference between previous research estimations of the groundwater influx variations and model’s results for the groundwater inflow into the lake, but the behavior prediction of the model is reasonable.

The data for domestic, agricultural and industrial water demand was derived from the Department of Water Resources and Irrigation at the Ministry of Agriculture and Land Reclamation of the Kyrgyz Republic. The initial population, immigration rate, emigration rate, birth rate, death rate, and population growth rate were collected for the Issyk-Kul Oblast. The types of land use/land cover change in the study area were classified as agricultural, natural and urban, and the data available for 1990, 2000, and 2010 was obtained from the Data-Center Sharing Infrastructure of Earth System Science, and the National Statistical Committee of the Kyrgyz Republic. Other important parameters, like the urban and rural population water demand per capita, net irrigation water demand per hectare, and water demand for industrial products were derived from the CAWATER Info-Indicators of Sustainable Development for Central Asia Countries (Kyrgyzstan water resources, Socioeconomic indices 1990–2010).

3.3. Model Calibration and Validation

Model calibration is a process by which certain model variables are adjusted to obtain a match between model output and historical data [66]. The model must be validated in terms of its structure, behavior, and in order to check its applicability and accuracy. First, a behavioral replication was used as a verification method to test whether the model can reproduce both qualitatively and quantitatively the behavior of key parameters [67]. The linear correlation coefficient (CF), Nash-Sutcliff efficiency coefficient (NSE), multiplicative bias (MBias), values of the mean absolute error (MAE), and root mean square error (RMSE) were calculated based on the difference between measured historical data and the simulated data [68,69]:

$$CF = \frac{{{\displaystyle \sum }}_{i=1}^n[(V_{\mathrm{s}} - \overline{V_{\mathrm{s}}})(V_0 - \overline{V_0})]}{\sqrt{{{\displaystyle \sum }}_{i=1}^n[{(V_{\mathrm{s}} -\overline{V_{\mathrm{s}}})}^2] {{\displaystyle \sum }}_{i=1}^n[{(V_0 - \overline{V_0})}^2]}}$$

(16)

$$NSE=1- \frac{{{\displaystyle \sum }}_{i=1}^n{\left(V_{\mathrm{s}} - V_0\right)}^2}{{{\displaystyle \sum }}_{i=1}^n{\left(V_0 - \overline{V_0}\right)}^2}$$

(17)

$$MBias = \frac{{{\displaystyle \sum }}_{i=1}^n{V}_{\mathrm{s}}}{{{\displaystyle \sum }}_{i=1}^n{V}_0},$$

(18)

$$MAE = \frac{1}{n}{\displaystyle \sum }_{i=1}^n{\left|\right(V}_{\mathrm{s}}{ - V}_0)|,$$

(19)

$$RMSE = \sqrt{\frac{1}{n}{\displaystyle \sum }_{i=1}^n{{(V}_{\mathrm{s}}{ - V}_0)}^2}$$

(20)

where V₀ = Observed value of the variable, V_s = Simulated value of the variable, $\overline{V_0}$ = Mean of the observed value of the variable, and $\overline{V_{\mathrm{s}}}$ = Mean of the simulated value of the variable.

The CF is used to assess the agreement between the simulated data and observed data. The range of CF values is between −1 to +1. A CF value of exactly +1 indicates a perfect positive fit, while a value of exactly −1 indicates a perfect negative fit. The values of NSE coefficient ranges roughly from −∞ to 1. A perfect estimation would result in an MBias value of 1. Underestimation will lead to values smaller than 1 and overestimation to values greater than 1. MAE and RSME are used to assigns the average magnitude of the error. The optimal values of the MAE and RSME are 0 [70].

In a second step, using a sensitivity analysis, the most important variables have to be identified for the developed model or system development. This means that the sensitivity analysis can identify the main links between observations, model inputs, and predictions. Thus, the variable selection is very important as they are the basis for the regulation strategy [71]. Calculating model outputs by changing one variable while keeping other parameters unchanged, the calculation formula is as follows:

$$S_q= |\frac{\mathrm{∆}q\left(t\right)}{q\left(t\right)}\times \frac{y\left(t\right)}{\mathrm{∆}y\left(t\right)}|,$$

(21)

where t is time, q(t) indicates the system’s state at time t, y(t) denotes the system parameters that affect the system’s state at time t, S_q is the sensitivity degree of state q towards parameter y while ∆q(t) and ∆y(t) represent increments of state q and parameter y at time t, respectively. The results obtained can help to identify the parameters that maximally affect the system’s behavior and thus they can potentially help guide the analysis of simulation scenarios.

4. Results

4.1. Water Level Variations of Issyk-Kul Lake

The annual average water level of Lake Issyk-Kul has fluctuated since 1980 (Figure 4), showing a downward trend throughout the 1980s and 1990s. The water level decreased from 1606.73 m in 1980 to 1606.18 m in 1998, dropping 0.57 m in 18 years (−3.17 cm/year), while from 1999 to 2011, the water level increased. The year 2012 was characterized by another small decrease with a water level of 1606.99 m at the end of that year (+0.81 m in 14 years; +5.79 cm/year). The highest water level during the study period was detected in 2011 (1607.02 m), and the lowest water level was recorded in 1998 (1606.18 m), while the years from 1998 up to 2011 registered an increasing lake level amplitude of 0.84 m, indicating a high temporal dynamic. These fluctuations of the Issyk-Kul water level are the result of the climate change and the anthropogenic water use in the catchment.

4.2. Factors Influencing the Water Balance

4.2.1. Climate-Hydrologic Change

According to the water balance equations, the lake water level of Lake Issyk-Kul is closely related to climatic-hydrological factors, including the precipitation within its basin, evaporation over the lake’s surface area, and surface water runoff. The annual precipitation, evaporation and natural runoff variability during the period of 1980–2012 are shown in Figure 4. The trends of surface runoff, precipitation and evaporation variations are basically consistent with the changing trend of the lake water level, and show the same pattern of decrease, increase, and decrease again. Over the period from 1980 to 2012, the mean annual natural runoff was 2421 million m³, the annual maximum and minimum surface runoff were 2888 million m³ and 1998 million m³, respectively. The mean annual precipitation during the study period was 313.46 mm. The annual maximum and minimum precipitation were 396 mm and 238.5 mm, respectively. And the mean annual evaporation from lake surface area was 933.96 mm.

4.2.2. Anthropogenic Activity

In the present study, human activities are characterized by the total population, land use changes and the gross industrial product. Kyrgyzstan’s economy is growing fast and an impact on the hydrology of the lake basin can be expected. The population of the Issyk-Kul Oblast as a proportion of the total population of Kyrgyzstan has declined slightly over the study period from 10.05% in 1980 to 8% in 2012, but the total population increased from 355,000 in 1980 to 448,000 in 2012. This equals an annual population growth rate of +0.712% for recent 33 years. This increase of the population could lead to an increasing use of water as the population change is the main driving factor behind the amount of water consumed in the lake basin. It not only determines the domestic water demand but also industrial water demand and the area of agricultural land, which in turn affects the agricultural water demand.

Land use patterns in the study area are divided into three categories, agricultural land, natural land, and urban land. The distribution of land use in 1990 was 9.9% agricultural, 88.7% natural, and 1.4% urban. In 2012, this changed to 10.5% agricultural, 87.9% natural, and 1.6% urban. Overall, the changes during this period were only small, with only minor changes in the water consumption. The analysis of the land use changes shown a process of ongoing urbanization in the Issyk-Kul catchment. Natural land and agricultural lands were found to be converted into urbanized areas by 0.012 and 0.009 ha per person. But the land use distribution alone stands for only one part of the anthropogenic influence on the water balance as the land use intensity is also a major factor. The annual water consumption for agricultural irrigation, domestic and industrial purposes has decreased from 1120 million m³, 16.95 million m³, and 10.59 million m³ in 1980 to 380 million m³, 6 million m³, and 6.8 million m³ in 2012, respectively. This is a drop in water consumption of 66% for the agrarian sector, 64% for the domestic sector and 36% for the industrial water use. This is the result of the disintegration of the Soviet Union between 1989 and 1992, which led to a decline of the agrarian and industrial water use, in combination with adopted lake water management measures. However, with the growth of the population in the lake basin, the local water demand has been rising, leading to a water level decline. This is being detected and controlled through various government policies, such as the use of other lakes and basins, management of lake streams, soil erosion prevention measures and rainwater collection [32,72]. Although the overall demand for agricultural water has been reduced, the water consumption of this sector is linked to the crop types grown in the basin and there is a great likelihood of an increasing water demand as more irrigation intensive crop types (e.g., cotton) are cultivated. In addition to the crop types, the applied irrigation techniques such as flood, furrow, sprinkler or drip irrigation have a significant impact on the water use efficiency and are related to processes like secondary salinization, which in turn contributes to an increase of the water demand due to leaching. Therefore, considering the best irrigation practices is conducive to the water management in the Issyk-Kul lake basin [73,74].

The GIP arises from industrial activity which requires water, though, in the case of the industry, much of the water is not “consumed” but recycled into the river network as wastewater. It could be argued that as economic development proceeds, such recycling will increase, as will the quality of the recycled water [75]. After Kyrgyzstan’s independence, the annual industrial effluent discharge fell from 674 million m³ in 1991 to 466 million m³ in 2000 (a decrease of 208 million m³ or 31%) due to the collapse of the industrial sector. The amount of treated industrial wastewater remained stable at approximately 150 million m³, which means that only 32% of the industrial wastewater was channeled through wastewater treatment plants. Other sewage discharges, mainly municipal wastewater and farmland drainage water, increased from 2.90 million m³ (1991) to 3.80 million m³ (2000). Only about 56% of Kyrgyzstan’s cities and towns are connected to a centralized sewage treatment system. More than half of the small and medium-sized cities, villages and towns have no centralized wastewater treatment, but are responsible for 27% of the total wastewater, which is discharged untreated.

4.3. Calibration and Validation

In this study, calibrations of the model provided new estimations of the groundwater inflow of the lake. The estimated coefficient of infiltration and water losses from groundwater using SELLA’s graphical function for calibration between 0.32 to 0.54, and 0.13 to 0.21, respectively. The model was validated using multiple tests, as described by Bayer [76], including structure assessment, dimensional consistency, extreme condition test, and behavior reproduction tests. The structure assessment test compares the model structure with the real system [55]. The dimensional consistency test ensures uniformity in the units of measurement of all the variables in the model [48]. An extreme condition test evaluates the model’s responses under extreme input values (see Section 4.3.1). The behavior reproduction test checks the model’s ability to replicate the behavior of the real system (see Section 4.3.2). The model performed satisfactorily under all these tests. The model was validated by comparing the model output with the observed data for the lake water level, the agricultural water demand, and the total population in the Issyk-Kul Oblast between 1980 and 2012.

4.3.1. Extreme Condition Test

The structure of a system dynamics model should permit extreme combinations of levels (state variables) being properly represented in the system. By examining the model structure for extreme conditions, the confidence in the model’s ability to behave plausibly for a wide range of conditions is developed and thereby the model’s usefulness to explore policies that move the system outside of historical ranges of behavior is enhanced [77]. The extreme condition used for this test in this study was to reduce the evaporation from the lake surface (the only outflow of the lake) to zero. As a result of the absence of evaporation the lake’s volume should rise progressively, which the model properly represented (Figure 5).

4.3.2. Behavior Reproduction Test

For the behavior reproduction test, the generated model behavior is judged against the recorded historical behavior. According to the Equations (16) to (20), the model was validated by comparing the simulation results with the historical data collected from 1980 to 2012. Based on the research objectives, the lake water level, agricultural water demand, and the total population in the Issyk-Kul Oblast were selected as significant variables for the validation of the SD model. The results of this model validation are shown in

Table 1. Analysis of the model validation for simulated and observed values of variables.

Variables	CF	MAE	RMSE	MBias	NSE
Lake level	0.928	0.0667	0.0906	1.000	0.9132
Total population	0.999	0.0002	0.0006	1.000	0.9998
Agricultural water demand	0.966	0.5816	0.7761	0.988	0.9327

and in Figure 6. They show the comparison of the simulation results and the observed data of the lake water level.

The validation results of the observed and simulated trends of the lake’s water level, agricultural water demand, and total population were found to be acceptable for a complex model. This suggests that the developed SD model is capable of reproducing the behavior of different parameters within the system and using this model is appropriate in this research.

4.3.3. Sensitivity Degree Analysis

In this study, to assess the degree of sensitivity for each state variable, selected parameters were increased or decreased by about ±10% for the study period from 1980 to 2012. Based on Equation (21), nine important model parameters were selected for a sensitivity analysis of the lake: precipitation (P); surface water flow (SW); groundwater inflow (GW); and evaporation (E) represented the natural factors; and population (p); agricultural water demand (A_wd); industrial water demand (I_wd); domestic water demand (D_wd), and; agricultural land area (L_agr) were used as anthropogenic factors. The water level of the lake was selected as the target variables for the sensitivity analysis.

Table 2. Sensitivity analysis of the model (%).

Parameter Change (%)	SW		GW		E		P
	+10	−10	+10	−10	+10	−10	+10	−10
Lake Level	1.47	1.09	0.57	0.63	−1.02	−0.53	0.41	−0.08
Parameter Change (%)	p		Awd		Iwd		Dwd		Lagr
	+10	−10	+10	−10	+10	−10	+10	−10	+10	−10
Lake Level	0.045	−0.040	0.16	0.08	0.006	0.007	0.007	−0.007	0.073	0.10

shows the results of this model analysis and the influence of variations in the selected parameters on the lake level. The variation rankings for each selected parameter are as follows: SW > E > GW > P > A_wd > L_agr > p > D_wd > I_wd.

Overall, the variation in the values of the agricultural area, population, domestic water demand, and industrial water demand had little impact on the lake water level. On the other hand, the water level was much more sensitive to the dynamics in the surface water flow, evaporation, groundwater inflow, precipitation, and the agricultural water demand. This attests that the inflow and outflow of water plays an important role for the water balance, and the impact of natural factor variability is still stronger than the impact of the human activities on the lake water level.

4.4. Simulation Scenarios and Results

The simulation includes a base run and four scenarios. The base run was used to estimate the lake water level from 2013 to 2040, and to assess the dynamics of water level vulnerability using the average values of the observation period between 1980 and 2012. For this simulation run, it was assumed that the averages of the observed values will remain unchanged in the future. Besides this base run, four scenarios were evaluated, each of them focusing on a different combination of the main driving forces for the lake water level variability (

Table 3. Different scenarios and description.

Scenario	Description
Base run	Variables used the average values of the observation period after 2012
River flow	Surface runoff assumed to close to maximum runoff and minimum runoff in past 33 years, respectively, after 2012
Water cycle components	Based the evaporation data in 2012, assumed the evaporation will increase by 1.76 mm per year after 2012; Precipitation changes under RCP 4.5 and RCP 8.5 scenarios after 2012
Decrease in water consumption	Assumed the water consumption decreased by 10% after 2012

). Figure 7 indicates the behavior of the water level variations under the extreme scenarios in the simulation period. The future scenarios for estimating the water level cover a period from 2013 to 2040; analyzing the relationships between the corresponding parameters are expressed quantitatively.

Scenario 1: Impact of the river runoff on the lake level:

The first scenario considers the impact of the river discharge into Lake Issyk-Kul based on historical data for 16 rivers. The historical discharge data series for the 16 main tributaries in the Issyk-Kul basin show that those rivers account for a total annual average runoff of 2421 million m³, with a maximum runoff of 2888 million m³ and minimum runoff of 1998 million m³ between 1980 to 2012 (Figure 4). The first scenario was run under the assumption that the river runoff in the future remains close to its maximum runoff and minimum runoff during the past 33 years.

Scenarios 2 and 3: Impact of the water cycle components on the lake level:

The historical annual average evaporation from the Lake Issyk-Kul surface from 1930 to 1950, 1951 to 1970, 1971 to 1990, and 1991 to 2010 was 820.31 mm, 884.81 mm, 907.05 mm, and 941.51mm, respectively. The evaporation increased from 800.34 mm in 1930 to 942.76 mm in 2012, increasing 142.42 mm in 83 years (+1.716 mm/year) [16]. The second scenario was run under the assumption that this positive trend continues for the period 2013 to 2040.

Scenario 3: The third scenario uses the annual evaporation increase from scenario 2 in combination with predicted precipitation changes. These predictions were based on the representative concentration pathway (RCP ) 4.5 and RCP 8.5 climate change scenarios and were extracted from the 20-km Weather Research and Forecasting (WRF) model (downscaled and bilinear interpolated) [78] for four meteorological stations (Balykchy, Karakol, Kyzyl-Suu and Cholpan-Ata) for the period 2013–2040.

Scenario 4: Impact of the water consumption on the lake level:

The annual water consumption for agricultural irrigation, domestic and industrial purposes has decreased from 1120 million m³, 16.95 million m³, and 10.59 million m³ in 1980 to 380 million m³ (−66%), 6 million m³ (−64%), and 6.8 million m³ (−36%) in 2012, respectively. This is the result of the disintegration of the Soviet Union between 1989 and 1992, which led to a decline of the agrarian and industrial water use, in combination with adopted lake water management measures which have increased the water use efficiency. This fourth scenario was run under the assumption that the water consumption will be decreased by 10% until 2040.

The simulated results from these four scenarios indicate that the water level fluctuations vary based on the applied water resource scenario. The scenario simulations (Figure 7) demonstrate that water cycle components and river inflow to the lake are the crucial factors. The behavior of the lake water level in response to changes in the model inputs indicate that:

• The behavior of the lake water level under the scenarios decrease in water consumption (scenarios 4) the annual lake water level will not change much; the annual average lake water level would increase from 1606.532 m to 1606.688 m (+0.156 m).
• The river flow and water cycle components scenarios (1, 2 and 3) have a greater influence on the water level. As the river runoff (scenario 1) close to its maximum runoff and minimum runoff during the past years, the annual average lake water level drops from 1606.532 m to 1605.460 m (−1.072 m), or increase to 1606.913 m (+0.381 m). Moreover, as the evaporation increases (scenario 2) by 1.76 mm during the future, the annual average lake water level drops from 1606.532 m to 1606.109 m (−0.423 m). Scenario 2, in combination with predicted precipitation changes (RCP 4.5 and RCP 8.5), as scenario 3, have a significant impact on the water level. This impact was greater than for scenario 2 and scenario 4, as the annual average lake water level would decrease from 1606.532 m to 1606.013 m (+0.519 m), or increase to 1607 m (+0.468 m).
• The analysis determines that the river inflow and the water cycle components scenarios significantly impact the lake water level. They are the key strategies for managing the lake water level. Reducing the water consumption in the catchment is indeed leading to an increase of the surface water inflow to the lake and in turn to a rising water level.

5. Discussion

The results presented here indicate that changes of the river discharge and the water cycle components are the main driving factors for the Issyk-Kul water level variations. The issue of water cycle components changes in response to the lake’s water level has been discussed for this region. Wang et al. [79,80] used the water balance equation and analyzed the factors of the water balance in the Issyk-Kul catchment and concluded that the correlation between these factors and the meteorological conditions is statistically significant. This assessment confirms the results provided by the SD model and presented here. Salamat et al. [16] studied data from hydro-meteorological stations around the Issyk-Kul and used correlation analyses to determine the fluctuations of the lake water level that are caused by surface runoff, precipitation, and evaporation. Burul et al. [44] collected data on water consumption to identify the main water consumers and their likely impact on Lake Issyk-Kul. The results of those two studies are based on linear correlation analyses and thus tend to provide less accurate insights than the results derived from the nonlinear SD approach used here. This is mainly because of SD model’s ability to provide proper insights into potential consequences of system perturbation, which are dependent on efficiently recognizing the main constituents and feedback loops among them. The study presented here takes the two main driving forces—water cycle component factors and anthropogenic factors—into account and considers the interactions between them and their influence on the Issyk-Kul water level. These non-linear correlations and quantitative analyses show that the water cycle components are an essential parameter for the dynamic of the lake water level.

Many researchers have developed SD models for lake water balance assessments in arid regions, and have demonstrated that the surface water runoff in arid regions varies both temporally and spatially and is the largest term in the lake’s water balance. Among various socio-hydrologic factors, surface water runoff is most commonly related to the lake water level. In a complex water balance system, the SD model is a great way to solve these complex problems as it provides a unique framework for integrating different physical and socioeconomic systems that are critical to the watershed management. Several studies in lake basins have reported how the different social, economic and environmental subsystems interact [52,81,82,83]. They developed SD models for the environmental planning and the management of lake basins, and also to simulate the impact of climate change and variability.

According to the Issyk-Kul SD model (Figure 1), some important points have been achieved. The model used in this study includes, as a novel approach, hydrologic and socioeconomic sub-models, implementing the previous research for water level fluctuation of the Issyk-Kul, and finalized versions of other simple models from other regions. The results show a complex response between variables, such as the connections between the population growth and the positive response of the agricultural area, as farming is the main economic activity in the Issyk-Kul basin. A continued increase of these factors may have a dramatic negative effect on the lake water level. On the other hand, according to the four simulation scenarios, certain changes, like decreases in the water consumption (scenario 4), would stabilize the Issyk-Kul water level. The scenarios showed that the lake system is highly vulnerable to water cycle components (scenarios 2 and 3) and how much inflow from its tributaries is required for a stable water balance (scenario 1).

Overall, the increase of the surface water inflow into the lake or the decrease of water consumption for irrigation alone cannot guarantee a sustainable development for the Issyk-Kul basin. In addition, maintaining a higher or lower lake level has an impact on the social and economic development of this region. Thus, in order to reach a socioeconomically sustainable development, the anthropogenic water consumption, the ecological water use and the proportion of the industrial water must be managed holistically, for example through the promotion of water-saving practices throughout society or the improvement of water use efficiency.

6. Conclusions

A dynamic water balance model of the Issyk-Kul lake basin was developed for the period from 1980 to 2012, and the relations between all relevant and available environmental and socioeconomic factors were formulated by mathematical equations. The test run with historical data validated this model (CF, NSE, and MBias > 0.91; RMSE and MAE < 0.78). In a second step, the main factors impacting the lake’s water level were evaluated via sensitivity analysis and various water resource scenarios. Results based on the sensitivity analysis indicate that the socio-hydrologic factors had different impacts on the lake water level change (SW > E > GW > P > A_wd > L_agr > p > D_wd > I_wd), with the main influence coming from the water inflow dynamic, namely, the increasing decreasing water withdrawal from lake tributaries. Land use changes, population increase, and water demand decrease were also important factors for the lake water level variations. Results of four scenario analyses demonstrated that water cycle components changes as the evaporation and the precipitation and the variations of river runoff into the lake are essential parameters for the dynamic of the lake water level. Overall, the water balance of the Issyk-Kul is highly complex, and any attempt for a sustainable water resource management needs to address a wide range of relevant parameters, from land use change and water consumption and population dynamics.