Vulnerability of the Caspian Sea shoreline to changes in hydrology and climate

The Caspian Sea is the largest inland water body of the world, located in Eurasia between Iran, Russia, Kazakhstan, Turkmenistan, and Azerbaijan (Iranian National Institute for Oceanography and Atmospheric Science Studies (INIOAS) 2020). The volume and area of the Caspian Sea are approximately 78 000 km³ and 371 000 km² (excluding Garabogazköl Bay area which is around 18 000 km² located east of the sea) (Kosarev et al 2009). It has a salinity of approximately $12\, \,{\rm gl^{-1}}$, about a third of the salinity of most seawater (Zonn et al 2010). The SWL of the Caspian Sea is approximately 28 m below mean sea level and mean depth of the sea is about 230 m (Amante and Eakins 2009). Based on the Köppen-Geiger climate classification, the Caspian Sea basin has diverse climates, such as cold and humid (climate type Df) in the Volga River basin, warm temperate humid (Cf) in the south and west, and arid (BW) and semi-arid (BS) in the southeast and northeast (Kottek et al 2006). In winter, 20 000–95 000 km² of the sea is frozen in north (Kouraev et al 2004, Ivkina et al 2017). Annual precipitation in surrounding rain gauges of the sea varies between 130 mm (Cheleken in Turkmenistan) to 1900 mm (Anzali in Iran). Volga River contributes over 80% of the total discharge (Arpe et al 2000, 2013, Ozyavas et al 2010, Roshan et al 2012, Chen et al 2017), with mean annual flow of 250 km³ (Iranian National Institute for Oceanography and Atmospheric Science Studies (INIOAS) 2020). The sea is divided into three regions: Northern, Middle and Southern Caspian, which are geometrically different (Amirahmadi 2000). The Northern region has mean depth around 5–6 m and represents less than 1% of the sea's volume. The average depth of the sea in Middle Caspian reaches to 190 m (Dumont et al 2004). The Southern Caspian is the deepest region, with sea depth greater than 1000 m (Lahijani et al 2019).

SWL data for the Caspian Sea (figure 2(a)) were produced based on information retrieved from the Hydroweb database (Cretaux et al 2011) and published observations from tide gauges (Kostianoy et al 2014). In the Caspian Sea, in-situ observations of SWL are available for period 1840–2000 (Kostianoy et al 2014), while the Hydroweb database provides SWL data obtained using satellite altimetry from 1992 to present (Cretaux et al 2011). A depth map of the Caspian Sea was taken from ETOPO1 (Amante and Eakins 2009) and annual water body maps of the sea were produced by Normalised Difference Water Index (NDWI) using Landsat images (appendix A). The area-volume-level relationships in the Caspian Sea were determined based on area and SWL data time series for the sea (appendix B).

Zoom In Zoom Out Reset image size

The cumulative capacity of reservoirs in the Caspian Sea basin was estimated based on the global FAO AQUASTAT database (FAO 2014). AQUASTAT gathers detailed information about dams in each country, especially on location, height, reservoir capacity, surface area and main purpose and it covers information on over 14000 dams. As shown in figure 2(c), most of dams in the basin were constructed in the 1950s and total capacity of reservoirs is 223 km³ which mostly are in the Volga River watershed. Comparing cumulative reservoir capacity and SWL revealed that, before the 1980s, there was a connection between decreasing SWL and increasing volume of reservoirs in the basin. Without river regulation, one study in the 1990s estimated that SWL in the Caspian Sea would be at least 1–1.5 m above the level observed at that time (Georgievsky and Shiklomanov 1994). Based on AQUASTAT (FAO 2014), the majority of dams on the Volga River and Ural River are single-purpose, for hydroelectricity, so oversupply flow must pass unused from such dams.

To estimate total inflow to the Caspian Sea, we used flow observations from downstream river gauges on major rivers taken from the Global Monthly River Discharge Data Set (RivDIS). RivDIS contains monthly averaged discharge measurements for 1018 stations located throughout the world (Vorosmarty et al 1998). The flow data included a considerable number of missing values after 1983, except for the Volga River at Vernelebyazhye gauge. Therefore, total inflow to the Caspian Sea was estimated based on the observed discharge at Vernelebyazhye in 1940–2015. Historical data obtained from RivDIS showed (appendix C and supplementary materials) that the Volga River supplies more than 80% of total inflow to the sea (Arpe et al 2000, Ozyavas et al 2010, Chen et al 2017).

Water body maps of the Caspian Sea were produced using NDWI (Gao 1996, Mcfeeters 1996):

$$\begin{equation}NDW{I_{NIR - SWIR}} = \frac{{\left( {NIR - SWIR} \right)}}{{\left( {NIR + SWIR} \right)}}\end{equation} \tag{ 1 }$$

$$\begin{equation}NDW{I_{Green - NIR}} = \frac{{\left( {Green - NIR} \right)}}{{\left( {Green + NIR} \right)}}\end{equation} \tag{ 2 }$$

where NIR, SWIR and Green are near infrared (wavelength 0.77–0.90 µm), shortwave infrared (wavelength 1.55–1.75 µm) and green (wavelength 0.52–0.60 µm) respectively in the electromagnetic spectrum. The main NDWI is based on Green and NIR bands showed in equation (2) (Mcfeeters 1996) and the rest equations are the modified version.

The Caspian Sea is in cloudy region, so we missed some satellite images. We used the Landsat images which had cloud cover $< 20{\text{% }}$ in studied time span from 1977–2018. In order to produce annual map of NDWI, we used summer images of each year when cloud cover is less compared to other seasons; then, for specifying water from land, NDWI values greater than zero were considered (Gao 1996, Mcfeeters 1996) as threshold. In appendix A, more details are presented on the procedure of producing NDWI maps using Google Earth Engine (Gorelick et al 2017) JavaScript Application Program Interface and the selection of NDWI threshold. We could not produce NDWI maps for 1978–1986 due to the low quality of images. We calculated NDWI for the period 1987–2018 using equation (1) based on the spectral resolution of recent Landsat sensors (ETM, ETM + and OLI) which cover $NIR$ and $SWIR$ bands, but we calculated NDWI of 1977 using equation (2) due to limited bands in the Landsat MSS collection spectral resolution which has only $Green$ and $NIR$.

Conceptually, the vulnerability of surrounding countries of the Caspian Sea is related to the importance of the shorelines in respect to the impacts of the SWL fluctuation. The vulnerability ratio, ${v_i}$ in equation (3) was considered to be a representative index for comparing the vulnerability of countries:

$$\begin{equation}{v_i} = \frac{{{A_i}}}{{{L_i}}}\end{equation} \tag{ 3 }$$

where ${A_i}$ is potential vulnerable area in each country, determined from comparing the NDWI maps of the sea in 1995 and 1977 when the sea had the highest and lowest areas respectively. In 1977, SWL (−29 m) was the lowest over the last 400 years (Kosarev et al 2009). ${L_i}$ is the length of coastline in each surrounding country. We used ${v_i}$ because the coastline of the Caspian Sea is not equally distributed, e.g. Kazakhstan has about 3000 km of coastline, while Azerbaijan has 650 km so vulnerability cannot be judged only by ${A_i}$.

We quantified the relationship between SWL raise and covered potential vulnerable area (CVA) by 33 NDWI maps from 1977 to 2018. According to 30-meter spatial resolution of Landsat, SWL-CVA relationship can be developed from small sites to a whole coastline running around surrounding countries. The SWL-CVA relationship was developed for the Caspian Sea, Northern/Middle/Southern Caspian and surrounding countries.

The annual SWL in the Caspian Sea was modelled based on the following equation:

$$\begin{align}{\text{ }}{L_{i + 1}} &= {L_i} + \frac{{Q_{i + 1}^v}}{{0.86 \times {A_{i + 1}}}} \times 1000 + \frac{{\left( {{P_{i + 1}} - {E_{i + 1}}} \right)}}{{1000}}\nonumber\\ &= {L_i} + \left( {Q_{i + 1}^t + {P_{i + 1}} - {E_{i + 1}}} \right)/1000\end{align} \tag{ 4 }$$

where ${L_i}$ is simulated SWL (m) in the Caspian Sea in year $i$ ($1979 \leqslant i \leqslant 2015$), ${A_{i + 1}}$ is the area of the sea in km² (excluding area of Garabogazköl), which was estimated by NDWI, ${P_{i + 1}}$ and ${E_{i + 1}}$ are precipitation (mm) and evaporation (mm) over the sea from the National Center for Environmental Prediction Climate Forecast System Reanalysis (NCEP-CFSR) model (Saha et al 2010). CFSR model available after 1980, which showed high accuracy when compared against in-situ data (appendix D). $Q_{i + 1}^v$ is Volga River flow in km³ at Vernelebyazhye in 1980–2015, and 0.86 is the proportional ratio of Volga River contribution to total flow to the sea. Total inflow is $Q_{i + 1}^t$ in mm shown in figure 2(b). Volga contribution is calculated by minimising the root mean square error (RMSE) between modelled and observed SWL and using that, the simulation results showed a good match between modelled and observed SWL ($correlation = 0.96$ and $RSME = 0.21$) (figure 2(a)). It should be noted that we considered Volga contribution equals to 0.86, in order to estimate the total sea inflow from all rivers . This approach was adopted due to the lack of available inflow data from other rivers (except Volga) and outflow from the sea (e.g. to Garabogazköl Bay).

The average inflow and net precipitation in different periods are shown in figures 2(b) and (d) by solid lines. In 1977, 1995 and 2008 (figure 2(a)), a significant change observed in the trend of level, inflow and net precipitation. After dam construction on the Volga (mainly in the 1950s), the flow of this river showed a declining trend from the early 1960s to the late 1970s (figure 2(b)). However, no dam has been constructed on the Volga River since the 1970s (figure 2(c)), and mean annual Volga inflow to the Caspian Sea has increased significantly since then. Mean net precipitation ($P - E$) from 1977 to 1995 increased from −755 to −692 mm. This increase was evident as a continuous raise in SWL between 1977 and 1995. From 1995 to 2010, net precipitation decreased, so SWL started to decline. After 2008, mean annual inflow has also declined, in parallel with continuing decreases in net precipitation, so SWL decline has become faster.

For global sensitivity analysis we applied method of Morris. Method of Morris (Morris 1991) is a so-called one-step-at-a-time method (OAT), meaning that in each run only one input parameter is given a new value that provides information regarding the overall effect of each input parameter on model output responses and the higher-order effects, such as interactions between parameters and non-linearity. For example, consider the model $M$ with a vector of $k$ parameters $\left( {{{{\theta }}_{i,j}} = 1, \ldots ,{\text{ }}k} \right)$ within the feasible parameter space ${{\theta }}$, that simulates $m$ response vectors of the system $\left( {{S_{i,j}} = 1, \ldots ,{\text{ }}m} \right)$:

$$\begin{equation}[{S_j}, \ldots ,{S_m}] = M\left( {{\theta _1}, \ldots ,{\theta _m}} \right).\end{equation} \tag{ 5 }$$

After running model $M$ for the given parameter sets, the local sensitivity measure (also referred to as the elementary effect, $EE$) is then computed for each parameter $i$ for model response $j$ as follows:

$$\begin{equation}E{E_{i,{\text{ }}j}}\left( \theta \right) = \frac{{{S_j}\left( {{\theta _1}, \ldots ,{\theta _{i - 1}},{\text{ }}{\theta _i} + \Delta , \ldots ,{\theta _k}} \right) - {S_j}\left( \theta \right)}}{\Delta } \end{equation} \tag{ 6 }$$

where $\Delta$ is a value in the predefined increments and $({{{\theta }}_{i,j}} = 1, \ldots ,{\text{ }}k)$ is a random sample in the parameter space so that the transformed point $\left( {{{{\theta }}_1}, \ldots ,{{{\theta }}_{i - 1}},{{{\theta }}_i} + \Delta , \ldots ,{{{\theta }}_k}} \right)$ is still within the parameter space θ. The resulting distribution $E{E_{i,{\text{ }}j}}{\text{ }}$ associated with each parameter ${{{\theta }}_i}$ is then analyzed to determine $\mu$, the mean of the distribution which assesses the overall importance of the parameter on the model output; and $\sigma$, the standard deviation of the distribution, which indicates non-linear effects and/or interactions. For non-monotonic models, some $EE$ values with opposite signs may cancel out when µ is calculated, and hence Campolongo et al (Campolongo, Francesca and Saltelli 1997) proposed the use of ${\mu ^*}$, the sample mean of distribution of absolute values of the elementary effects. Here we introduced evaporation, precipitation and total inflow as parameters (θ) of SWL simulation model ($M$) and average of SWL in all years of simulation as the response of the model ($S$).

According to Torabi Haghighi et al (2016), closed lakes can reach an equilibrium state in response to given hydro-climatological conditions. This means that if the lake water balance components in any period (e.g. 20 years) repeat infinitely, after a response lag time, the fluctuation in lake water level, volume and area can converge to an equilibrium state. However, in some situations, level can keep increasing (e.g. forming open lakes) or decreasing trend (e.g. desiccated lakes). In the equilibrium state, lakes reach a specific area where the volume of inflow becomes equal to the volume of negative net precipitation (${Q^t} + P - E = 0$) so SWL becomes constant (Szesztay 1974, Mason et al 1994, Crétaux and Birkett 2006, Haghighi and Kløve 2015, Torabi Haghighi and Kløve 2017, Torabi Haghighi et al 2018).

To assess the equilibrium state of the Caspian Sea, we defined two sets of change factors, ${\alpha _Q}$ and ${\alpha _{P - E}}$ (ranging between 0.9 and 1.1 with 0.01 increments), for inflow and net precipitation. Overall, a total of 421 ($21 \times 21$) possible scenarios were generated to consider combinations of ${\alpha _Q}$ and ${\alpha _{P - E}}$. For each scenario, time series from 1980 to 2015 of annual $P - E$ (from CFSR) and ${Q^t}$ (total inflow) were multiplied by ${\alpha _{P - E}}$ and ${\alpha _Q}$ respectively. In each scenario, multiplied $P - E$ and ${Q^t}$ were repeated infinitely; then, SWL was simulated. In some scenarios, simulated SWL converges to constant value which is equilibrium SWL. Finally, we plotted equilibrium SWL heatmap in which each combination of probable inflow and net precipitation can be mapped to a specific equilibrium SWL. Using the SWL-CVA relationships and equilibrium SWL heatmap, we created CVA heatmaps in equilibrium states for the Caspian Sea, Northern/Middle/Southern Caspian, as well as surrounding countries.

We also determined the sensitivity of equilibrium SWL to changes in inflow and net precipitation. For this purpose, we fixed net precipitation and investigated the relationship between inflow and equilibrium SWL, and then fixed inflow and investigated the relationship between net precipitation and equilibrium SWL. The slope of these relationships shows the sensitivity of equilibrium SWL to changes in inflow or net precipitation.

Based on NDWI maps, the minimum area of the Caspian Sea over the last 400 years was estimated to be 355 000 (in 1977) and the maximum area was estimated to be 380 000 km² (in 1995) during the past century (figure 3(a)). Therefore, the area of the Caspian Sea potentially vulnerable to desiccation was calculated to be 25 000 km². In 1977–1995, SWL has increased about 3 m and determined potential vulnerable area is consistent with depth less than 3 m (figures 4(a) and (b)). Although this area is not the absolute historical difference between maximum and minimum area in the Caspian Sea, it reflects the most recent (from 1940 to present) observed variation in this water body. Around 70% of the potential vulnerable area is located in the Kazakhstan territory, followed by Russia (5055 km²), Iran (1089 km²), Azerbaijan (578 km²) and Turkmenistan (482 km²) (figure 3(b)). The most severe coastline retreats occur in the Northern Caspian region bordering Kazakhstan and Russia (figure 3(a)). In Kazakhstan, the coastline retreated by more than 160 km in 1977 compared to 1995. In addition, for each 1 km of coast in Kazakhstan, more than 6 km² of coastal area are potentially vulnerable. The vulnerability ratio, $\nu$ in equation (3), for whole Caspian Sea and for its coast in Russia, Iran, Azerbaijan and Turkmenistan is around 4, 4, 1.5, 1 and 0.5 ${\text{k}}{{\text{m}}^2}$ km^–1 respectively (figure 3(c)).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Three ecologically important bays on the eastern coast, Dead Kultuk, Türkmenbaşy and Miankale, are highly vulnerable to SWL fluctuation (figure 3(a)). As shown by purple curve in figure 3(a), the majority of potential coastline retreatment in Kazakhstan is in Dead Kultuk and potential vulnerable areas in Iran are mainly located in the Miankale Gulf. Marginal shallow gulfs in the south of Russia, centre of Turkmenistan and south of Azerbaijan are also vulnerable to coastline retreat compared with past historical observations. For example, in Türkmenbaşy Bay, the sea retreat can be more than 25 km. In Azerbaijan, although the Bay of Baku shoreline has experienced no change during the past 80 years, the coastline retreatment in Ghizil-Agaj State Reserve is considerable. In the west of the Caspian Sea, the Kizlyar region is the most vulnerable area. Astrakhan Nature Reserve on the western shoreline is another highly vulnerable region. Around 90% of the potential vulnerable area of the sea is on the Northern Caspian coast, 10% is on Southern Caspian coast, and less than 100 km² (0.4%) is on the Middle Caspian coast.

Based on the available images from 1977 to 2018 (33 NDWI maps) and SWL data for this period, the relationship between SWL raise in m (MSL) and CVA in km² was calculated at different scales from the whole of the Caspian Sea to coastal regions in surrounding countries (table 1). The slope of the linear regression between SWL and CVA for the whole of the sea was $7684{\text{ k}}{{\text{m}}^2}$ m^–1 (figure 4(c)), which indicates that with a 1 m rise in SWL, about 7700 km² of the potential vulnerable area of the Caspian Sea would be covered. Linear approximations in table 1 are restricted to a period from 1977 to 2018 ($- 29 < SWL < - 26$) so, based on assumptions in regression models, any extrapolation of a fitted regression equation beyond the range of SWL can lead to biased estimates.

Table 1. Relationship between sea water level (SWL, m) and covered potential vulnerable area (CVA, ${\text{km}}$²) for the Caspian Sea and surrounding countries.

Region	Equation	r2
Caspian Sea	$CVA = 7684 \times SWL + 225781$	0.91
Northern Caspian	$CVA = 7103 \times SWL + 208261$	0.90
Middle Caspian	$CVA = 17 \times SWL + 570$	0.68
Southern Caspian	$CVA = 351 \times SWL + 11249$	0.56
Kazakhstan	$CVA = 6046 \times SWL + 175349$	0.88
Russia	$CVA = 352 \times SWL + 14141$	0.42
Iran	$CVA = 150 \times SWL + 4991$	0.65
Azerbaijan	$CVA = 36 \times SWL + 1519$	0.75
Turkmenistan	$CVA = 173 \times SWL + 4985$	0.66

Results of sensitivity analysis (figure 4(d)) indicate that the overall importance of evaporation (represented by ${\mu ^{\text{*}}}$) on the SWL simulation model output is highest followed by total inflow and precipitation respectively. In terms of $\sigma$, total inflow has the highest value followed by evaporation and precipitation that means total inflow is the variable involved in interaction with other variables and its effect is non-linear in SWL simulation model.

We found that each 1 km³ inflow alteration changes the equilibrium SWL by about 16–18 cm, while each 10 mm variation in P-E changes the equilibrium SWL by about 6–6.5 cm. It should be noted that the highest SWL in the Caspian Sea from 1980 to 2018, the period for which the SWL-CVA relationship is developed, equals to −26 m (in 1995). However, in the 1880s, $SWL = - 25.5$ m has been observed. Therefore, heatmaps of CVA for different parts of the sea (figure 5) were produced by extrapolating the SWL-CVA relationship for the range from −26 to −25.5. In this range, some coastal areas which were dry in 1995 will go under water and CVA can reach 27 500 km² in figure 5(b), i.e. more than potential vulnerable area (25 000 km²).

Zoom In Zoom Out Reset image size

The highest SWL in the Caspian Sea observed from 1840 to present was −25.5 m, so we considered equilibrium SWL more than −25.5 as one class (blue zone: full state, figure 5(a). Available records indicate that SWL less than −29 m has not occurred in the past, so equilibrium SWLs below −29 m are also categorised into one class in figure 5(a) (red zone: worst state). In 1987, when SWL was −27.7 m, the Dead Kultuk Bay was completely desiccated, shown in figure 3(a). We show this SWL in figure 5(a) as a white zone in the heatmap, before restoration work on Dead Kultuk Bay as an important environmental icon of the Caspian Sea. Therefore, around −27.7 m is a transient SWL between warning and restoring zones of the equilibrium SWL heatmap in figure 5(a).

Using calculated equilibrium $SWL$ (figure 5(a)) and $CVA - SWL$ regression model (table 1), we quantified CVA and water balance parameters (net precipitation and inflow) relation separately for surrounding countries of the sea in figure 6 As shown in this figure, Kazakhstan is the most vulnerable country by almost 18 000 km² of potential vulnerable area. In recent years from 2010 to 2015, the average of $P - E = - 786$ mm and ${Q^t} = 260$ km³. Under current situation, to recover all vulnerable area of Kazakhstan (reaching the contour of 18 000 km²), the total inflow should increase from 260 to more than 300 km³ or net precipitation should increase from −786 to −680 mm. Also, contour 18 000 km² determines other possible combination of $P - E$ and ${Q^t}$ for recovering all potential vulnerable area in Kazakhstan. As described before, to go beyond 18 000 km² and reach 20 000 km² contour, the SWL of the sea should increase from −26 m (highest observed SWL in 1995) to −25.5 m.

Zoom In Zoom Out Reset image size