**Shoreline Dynamics and Evaluation of Cultural Heritage Sites on the Shores of Large Reservoirs: Kuibyshev Reservoir, Russian Federation**

#### **Ionut Cristi Nicu 1,\*, Bulat Usmanov 2, Iskander Gainullin <sup>3</sup> and Madina Galimova <sup>3</sup>**


Received: 15 February 2019; Accepted: 20 March 2019; Published: 21 March 2019

**Abstract:** Over the last decades, the number of artificial reservoirs around the world has considerably increased. This leads to the formation of new shorelines, which are highly dynamic regarding erosion and deposition processes. The present work aims to assess the direct human action along the largest reservoir in Europe—Kuibyshev (Russian Federation) and to analyse threatened cultural heritage sites from the coastal area, with the help of historical maps, UAV (unmanned aerial vehicle), and topographic surveys. This approach is a necessity, due to the oscillating water level, local change of climate, and to the continuous increasing of natural hazards (in this case coastal erosion) all over the world. Many studies are approaching coastal areas of the seas and oceans, yet there are fewer studies regarding the inland coastal areas of large artificial reservoirs. Out of the total number of 1289 cultural heritage sites around the Kuibyshev reservoir, only 90 sites are not affected by the dam building; the rest had completely disappeared under the reservoir's water. The scenario of increasing and decreasing water level within the reservoir has shown the fact that there must be water oscillations greater than ±1 m in order to affect the cultural heritage sites. The results show that the coastal area is highly dynamic and that the complete destruction of the last remaining Palaeolithic site (Beganchik) from the shoreline of Kuibyshev reservoir is imminent, and immediate mitigation measures must be undertaken.

**Keywords:** cultural heritage; shoreline dynamics; GIS; UAV; Palaeolithic; Volga; European Russia

#### **1. Introduction**

The construction of large reservoirs along the large rivers of the world has, eventually, different effects: Local micro-climate modifications, disruption on the river flow regime [1], sediment transport [2,3], fauna [4], water chemistry [5], shore morphology [6–8], archaeology [9], fish yields [10], among other issues. They can also act as a place where different types of pollutants accumulate, and, in this way, it is easier to assess historical pollution [11]. One of the main effects is the triggering and the fast mechanic action of waves. These effects are accentuated by the global climatic changes, which are exponentially increasing every year.

Many studies deal with risk assessment [12,13], management [14,15], vulnerability [16–18], conservation strategies [19,20] and sustainability issues [21] regarding the cultural heritage of the coastal areas of seas and oceans. However, there is a lack of studies dealing with inland shorelines of large man-made reservoirs [22]. The Volga River is the largest river in Europe with a basin area of

1,360,000 km2; it is considered the main river in Russia, and its basin represents the most significant economic region in Russia [23]. During the Soviet Union, there was a usual practice to flood large territories in order to obtain electricity and to relocate a large number of inhabitants and their houses. Unfortunately, cultural heritage sites do not enter this category; they cannot be relocated or moved.

During the Soviet period (the late 1930s), the "Great Volga Scheme" was initiated; the purpose was the construction of a chain of dams along the Volga River and one of its major tributaries—the Kama River. The reservoirs of the Volga-Kama cascade are one of the largest cascades in the world, totaling 11 reservoirs (Figure 1, Table 1). The main purpose of the dams was to produce electricity; before the 1930s, the Volga was used only for transport and fishing [24,25]. As shown in Table 1, Kuibyshev reservoir has the largest surface and the highest number of types of uses.

There have been limited studies referring to the destruction of archaeological sites around the Kuibyshev reservoir [26,27], but there are no studies referring to the entire surface of the reservoir. Therefore, this study is necessary to assess the exact number of sites impacted by the reservoir creation in 1957 and to draw attention for local authorities in their mission for future management plans [28] of the shoreline area [29]. A detailed case study was chosen to demonstrate the destructive potential of wave erosion; this was accomplished by a systematic monitoring process. The main scope of this article is (1) to track the major changes of the Volga River after the construction of the Kuibyshev reservoir with the help of GIS (2) to identify the area(s) that contain the highest concentration of archaeological sites (3) to analyse how many archaeological sites were impacted following the construction of the reservoir (4) to monitor the evolution of the only left Palaeolithic site—Beganchik from the shores of Kuibyshev reservoir, which has been specifically chosen because of its high erosion rates.

**Figure 1.** Detail of the reservoirs of the Volga-Kama cascade.


**Table 1.** The main characteristics of reservoirs from the Volga-Kama cascade [24,25] (the numbers in the first column correspond to the reservoirs from Figure 1).

Legend: YOC—year of commissioning; RA—reservoir area; IC—installed capacity; AO—annual output; TOU—type of use; Fi—fishery, Fl—flood control, I—irrigation, Navigation, P—power production, R—recreation, T—timber rafting, W—water supply, Wr—water releases (sanitary, irrigation).

#### **2. Study Area**

Kuibyshev reservoir is a result of the construction of the Zhiguli Hydroelectric Station, Samara region, located between Zhigulevsk city (right bank of the Volga) and Tolyatti (left bank of the Volga); the reservoir covers the territory of regions Chuvash, Tatarstan, Ulyanovsk, and Samara. Kuibyshev reservoir has a surface of 6450 km2, a volume of water of 58 km3, a length of approximately 510 km, a mean depth of 9.3 m; these impressive numbers make it the largest reservoir in Europe [30], with a sedimentation rate of 8 mm/year. Important changes occurred in what concerns the sedimentation rate, which has fallen to 2.7–2.9 mm/year, after the commissioning of the dam in 1957. One of the main sources of sediments is a result of the abrasion processes, collapses of a huge amount of sediments into the reservoir [25,31].

Our area of interest is located in Tatarstan region (Figure 2a), on the left bank of Kuibyshev reservoir (Figure 2b), at the junction of Kama River in the Volga, about 75 km south-east of the city of Kazan (the capital city of Tatarstan). Beganchik site is located at approximately 2.8 km North-east of Izmeri village and 1.5 km north-west of Komintern village, on an isolated hill of the terrace above the floodplain; on the left bank of the confluence of Kama and Volga rivers, at the mouth of Aktai river (Figure 2c, Figure 3a).

The geology of the area consists of Permian, Pliocene and Quaternary deposits. Quaternary sediments are dominant in Volga-Kama terraces, eight palaeohydrological phases were identified from high and low fluvial activity; the most recent active phase corresponds to the Little Ice Age [32]. There is a limited number of studies regarding the evolution of the coastal area in the Tatarstan region, Russia [25], along with the analysis of landslides [33] and gully erosion [34].

**Figure 2.** Geographical location of: (**a**) Kuibyshev reservoir in a regional context; (**b**) Kuibyshev reservoir and the cultural heritage sites around it; (**c**) Beganchik site on the topographic maps scale 1:50.000 (edition 1984).

**Figure 3.** (**a**) Aerial image of Beganchik site from 2018; (**b**) Overlapping Volga River before and after the building of Kuibyshev reservoir, along with the inventory of archaeological sites.

#### **3. Archaeological Background**

#### *3.1. General Overview*

Ever since the Palaeolithic, large rivers and their fertile plains have been a magnet for prehistoric people to place their settlements; water is undoubtedly the most important resource that a community of people needs in order to decide where to place a settlement [35]. According to the local geographical factors and paleogeographic evolution of the landscape, the old location could be used by the next population of a different historic period. This is how multi-stratified archaeological sites were created. Around Kuibyshev reservoir, 1289 cultural heritage sites have been identified. The protection of cultural heritage assets in Russia was and is still a current issue, because of the country's complex political background; at present, cultural heritage is protected by the Federal Law 73-F3, On the Objects of Cultural Heritage (Historical and Cultural Sites) of the Peoples of the Russian Federation [36].

In this area, the oldest traces are attributed to Upper Palaeolithic-Mesolithic period. The area surrounding the Volga River has tremendous potential in regard to cultural heritage sites of Palaeolithic [37,38], Mesolithic age [39], Neolithic [40], Chalcolithic/Bronze Age [41], Early Iron Age [42], Middle Ages, etc. The only remaining Palaeolithic site which has not been impacted by the reservoir is Beganchik. Besides the Palaeolithic site, Beganchik, around the Kuibyshev reservoir there are many archaeological sites of international and national significance; among them, the Bolgar archaeological site. The Bolgar Historical and Archaeological Complex are part of the UNESCO World Heritage List since 2014; it represents the existence of the Volga-Bolgar civilisation (7–15th centuries AD), and the first capital of the Golden Horde in the 13th century [25].

#### *3.2. Beganchik Site*

Beganchik site was studied for the first time in September 1985 by M. Sh. Galimova and K. E. Istomin at the recommendation of E.P. Kazakov; the first description of the site is the islet named "the Izmeri Island". In 1981, mammoth fauna fossils were found; they were located on the towing-path in the south-western part of the islet, at the foot of the narrow, long butte, which had the shape of a peninsula with a length of about 200 m. The discovered mammoth fossils were five teeth and leg bones, together with large flint nuclei and tools in an area of 20 × 20 m2. Unfortunately, by the year 2000, this peninsula was completely eroded by the Kuibyshev reservoir. In the next years, almost every autumn (until 2012) Kazakov collected stone artefacts and faunal remains at the south-western tip of the islet in conditions of low reservoir level [43].

M. Sh. Galimova in 1985–1987 and 2000 also conducted investigations of this site. In 1986, a reference point was installed at the highest point of the islet, therefore all excavations and trenches were referenced after it. An excavation area of 104 m2 was located on the edge of the steep west coast of the islet; the cultural layer has been found at 100–130 cm depth; 1968 of artefacts were found. The specific features of the Beganchik stone industry, which was based on blade production by means of striking technique and its flint inventory, allowed M. Sh. Galimova to frame the site of initial (Upper Palaeolithic) period of the Mesolithic Ust-Kama culture. The main diagnostic tool of the Ust-Kama inventory is the arrowhead in a trapezoid shape with concave sides, which were shaped by retouching [44]. In 2000, the rescue excavations of the site were continued by M. Sh. Galimova with the participation of I. I. Gainullin. By that moment excavation territory of 1986–1987 was eroded by the reservoir. In general, the western and northern coasts of the islet were washed away by 20–25 m from the erosion ledge for 14 years (1986–2000). In the autumn of 2012, rescue investigation of the Beganchik site and Izmeri I site was conducted by the expedition of the National Centre for Archaeological Research of the Tatarstan Academy of Sciences. In the autumn of 2013, rescue investigation on the Beganchik islet was continued by a joint campaign of the "Expedition for Prehistory" of the Institute of Archaeology of the Tatarstan Academy of Sciences and "Archaeological Expedition" of the Chuvash State Institute for Humanities. The total excavated area was 20 m2; following the excavation and

surface findings, a significant collection of stone artefacts (439 items) and 80 bones of a mammoth were found [45].

#### **4. Materials and Methods**

In order to determine the shoreline dynamics, topographic map scale 1:300,000 (edition 1945), and Google Earth images from 2010, were employed. From the topographic map, the extent of the Volga River before the reservoir construction was digitised; from Google Earth images from 2010, the extent of the reservoir was digitised. They were overlapped in ArcGIS and the highest differences were observed.

The archaeological inventory, as a point feature, (Figure 3b) was provided by the Institute of Archaeology of Tatarstan Academy of Sciences. The database was compiled over a long period of time, both prior and after the filling of the reservoir; first sites were described in the early 1940s until the early 1960s, when special survey expeditions were undertaken, with the aim to find as many sites and record brief information about them. After the filling of the reservoir, more expeditions were undertaken to highlight the impact on the sites. Other sources in building the database included the descriptions of the Archaeological Maps of Tatarstan, Ulyanovsk and Samara regions. A survey has been unsystematic across the study area, sites are located to varying degrees of accuracy and the full extent of individual sites is not necessarily known. The database is still under construction, as the area is very large and only a few people within the Institute of Archaeology of Tatarstan Academy of Sciences is working to continuously update it. However, it is the most comprehensive archaeological dataset which currently exists in this area, hence will be used for this study. Density analysis of the settlements for the main historical periods was performed; this was made using the Point Density feature (using the circle as neighbourhood option) from ArcToolbox (ArcGIS).

The danger towards increasing and decreasing water level over the digital elevation model (DEM) was evaluated by making four working scenarios; the DEM used in this study is based on the Shuttle Radar Topography Mission (SRTM), with a pixel size of 30 × 30 m2. First, the water level was decreased by 0.5 m and 1 m, followed by increasing the water level with 0.5 m and 1 m, respectively. Changes in reservoir level regime occur out of two main reasons: Natural seasonal changes in the flow and artificial regulations of water discharge through hydraulic structures, the difference in baric pressure, wind speed and changes in the hydraulic slope. The water level of the reservoir is controlled by a special department—RusHydro. We chose these values taking into consideration our large-scale study area and to point out the minor oscillations in water level by arbitrarily increasing/decreasing it by ±0.5 to 1 m. In order to have a better image of the changes occurred along the Volga River after the Kuibyshev reservoir was built, the entire area was divided into three sectors, as follows: Sector 1 (Figure 4b), Sector 2 (Figure 4c), and Sector 3 (Figure 4d).

For monitoring the Beganchik site, a rich cartographic background, including Soviet aerial images from 1958 and 1980, Roscosmos aerial images (2008), and UAV flights [46] from the summer of 2017 and 2018 were used. As shown in numerous studies, old maps and aerial images are an important source of information [47] that can be easily digitised, integrated into GIS [48], and used in the field of cultural heritage [49,50]. All the maps and aerial photos have been georeferenced with the help of ArcGIS. The ground control points used for the drone flights were measured with a GNSS (Global Navigation Satellite System) receiver Trimble Geoexplorer 6000 XH and Leica Zeno 20. The UAV is a drone, model DJI Phantom 4. The survey was performed with a 12-megapixel camera mounted on the quadcopter; the UAV was controlled from a smartphone using Pix4D Capture software, which allows configuring the shooting parameters. Aerial photography was performed with the following parameters, height—70 m, picture overlapping—60–80%, camera position—90 degrees, meteorological conditions—no precipitation, and wind no more than 15 m/s. The photos were processed using the algorithms built into the Agisoft Photoscan software; the resulting model was processed by a polynomial approximation exponential kernel (PAEK) method with 1 m tolerance.

**Figure 4.** (**a**) General division of the Kuibyshev reservoir; (**b**) Sector 1; (**c**) Sector 2; (**d**) Sector 3.

#### **5. Results**

Each sector will be analysed according to the GIS integration of the spatial data collected from old maps and modern aerial images, followed by the analysis of the archaeological sites patterns and dynamics along the Volga River; then, the changes of Volga River will be analysed in the context of how many cultural heritage sites are directly affected by the reservoir construction. The four working scenarios will be analysed in order to evaluate the endangered sites towards increasing/decreasing water level of the reservoir. Finally, the monitoring results of the only left Palaeolithic site—Beganchik will be presented; this site has been specifically chosen because of its high erosion rates and being the only remaining Palaeolithic site around Kuibyshev reservoir.

#### *5.1. Volga River Dynamics*

From the town of Tver to Volgograd, Volga River flow velocity is affected by the 8 reservoirs. The reservoirs were built to control seasonal changes in flow; however, there are no significant changes when it comes to river discharge and the total annual discharge. In the middle Volga, the mean annual flow from 1876 to 1940 was 2876 m3/s; after the construction of reservoirs, from 1942–1955, the mean annual flow was 2780 m3/s [51].

**Sector 1** (Figure 5a) stretches approximately in the north-western part, next to the Zvenigovo city till south-east, at the junction between Volga and Kama rivers; it has a length of approximately 145 km. The most important cities within Sector 1 are Kazan (with a population of about 1,2 mil. people) and Zelenodolsk (with a population of about 98,000 people). Out of the three sectors, Sector 1 has the fewest changes, compared with the others; the most significant changes are located about 43 km downstream from Kazan city. Initially, Volga had a width of 1.4 km, while after the building of the Kuibyshev reservoir the width of Volga reached 9.5 km. Another significant change is located between Zelenodolsk and Kazan, at the junction of Sviyaga River in the Volga; from a width of 0.6 km, Volga reached a width of 11.2 km; except this, the reservoir water has mainly covered the left side of

the river. This is due to the geomorphological characteristics of the area; the right side represents the Volga uplands, while on the left side are the terraced plains of the lowland Volga River region [52].

**Sector 2** (Figure 5b) stretches from east to west on a distance of approximately 150 km and represents the lower Kama River junction to Volga River; before the Kuibyshev reservoir, the area located around the junction had a significant number of villages (which were completely destroyed). Moreover, important landscape changes occurred, along with an acceleration of coastal erosion with a direct effect on cultural heritage [25]. The most important city in this sector is Chistopol (with a population of approximately 60,000 people). Along its approximately 150 km length, after the building of Kuibyshev reservoir, Sector 2 had more or less a balanced development of the right and left bank; this is because both of the sides are located within the terraced plains of the lowland Volga River region. From an average width of 0.8–1 km, the Kama River reached widths of 13–36.8 km. From these numbers, we can realise the real proportions of the consequences of the Kuibyshev reservoir being built.

**Figure 5.** Results of the coastal dynamics analysis for: (**a**) Sector 1; (**b**) Sector 2; (**c**) Sector 3.

**Sector 3** (Figure 5c), with a length of approximately 263 km, stretches from Kama and Volga junction until the dam of the Zhiguli Hydroelectric Station, located between the cities of Zhigulyovsk and Tolyatti. The most important cities in this sector are Tolyatti (with a population of approximately 720,000 people), Ulyanovsk (with a population of approximately 614,000 people) and Bolgar (with a population of approximately 9000 people). Bolgar is well known for the Bolgar Historical and Archaeological Complex World Heritage site. Similar in development with Sector 1, the left side of the river being more developed than the right one; again, this is to the lower altitudes of the terraced plains of the lowland Volga River region [52]. Within this sector, the width of the Volga River before Kuibyshev reservoir was ranging from 0.7–2.1 km, and from 9.1–32.2 km after the building of Kuibyshev reservoir. Having this enormous width, it is sometimes called the Kuibyshev Sea.

#### *5.2. Archaeological Site Analysis*

Following the analysis of the archaeological database provided by the Institute of Archaeology of Tatarstan Academy of Sciences, the following periods were identified: Palaeolithic/Mesolithic, Neolithic, Chalcolithic/Bronze Age, Early Iron Age, Migration Period, and Middle Ages. Large river systems, e.g., the Volga, act as a magnet when it comes to taking a decision to place a prehistoric settlement. That is why, in the close proximity of the Volga River and its tributaries, there is a high density of archaeological sites. Water represents the main resource in establishing the placement of a settlement; this is documented and well-known across the archaeologists and geo-archaeologists [53].

On the basis of the existing database, the areas with the highest concentration of archaeological settlements attributed to a certain period will be identified and highlighted accordingly. Usually, the dynamics of the settlements are influenced by different factors, like climate change [54], natural

hazards [55] and threats from other populations. Having knowledge of the spatiotemporal distribution patterns of archaeological sites is a powerful tool to understand past human-environment interactions and to evaluate landscape vulnerability to natural [56] and anthropogenic changes [49].

As can be seen in Figure 6, the highest concentration of settlements for all the periods is located at the junction of the Kama and Volga Rivers, this is representing an important communication route. During the Palaeolithic/Mesolithic period (Figure 6a), the hunters-fishers-gatherers population was well adapted to the living around water bodies and forests; the highest concentration was at the junction of Kama and Volga rivers, followed by the adjacent areas of upstream and downstream of the junction. A good concentration can be observed on the Volga River, around the area where presently the city of Kazan is located; the thrive of settlements was due to the optimum climatic conditions for the Preboreal period [37,39], along with the highest levels of rivers and lakes, which was typical for the Mesolithic epoch [57]. The settlements were not very homogenous during this period. However, this can be observed during the Neolithic period (Figure 6b), when more settlements appear in the area between the today Kazan city and Kama-Volga junction. The Neolithic period is characterised by the emergence of pottery, new types of stone tools and the transition to sedentism with the help of active fishing and hunting. The majority of Neolithic sites are located on the remnants of the floodplain of the small rivers of the Kama River tributaries or on the first terrace of the Kama River [40].

**Figure 6.** Location and density analysis for archaeological sites for the following periods: (**a**) Palaeolithic/Mesolithic; (**b**) Neolithic; (**c**) Chalcolithic/Bronze Age; (**d**) Early Iron Age; (**e**) Migration Period; (**f**) Middle Ages.

Following Chalcolithic/Bronze Age period (Figure 6c), it can be observed even a higher degree of homogeneity among the settlements; this is due to the fact that the lowest levels of water were recorded in the Bronze Age. As a consequence of this, even the lowest altitudes were chosen to place the settlements, which is why the analysis shows a larger continuous surface

Figure 6d illustrates the density of the Early Iron Age settlements, which started to be more fragmented. The highest concentration is at the confluence of Kama-Volga Rivers and on the territory of today Bolgar, followed by scattered low-density areas the upstream Volga, at the mouth of Sviyaga River, and the downstream Volga. The settlements appear scattered because of their higher altitudinal position (higher position throughout the Holocene), due to the associated high flood levels [57]. As can be seen in Figure 6e, the Migration period is characterised by spreading of population downstream Volga River, until today Ulyanovsk city. However, the highest concentration is still located at the Kama-Volga junction; the fact that during this period the area is very poorly populated is also indicated by [58]. Finally, the Middle Ages (Figure 6f) show the highest fragmentation of the settlements. The highest concentration remains the same (Kama-Volga junction), while the settlements are scattered downstream and the upstream Volga until today Tolyatti and Zelenodolsk, respectively. During this period, the settlements are so scattered, due to the fact that the climatic conditions were suitable for the long-term occupation of river and lake floodplains. This is the turning point when people start to settle and make semi-permanent settlements and start off using the floodplain in order to practice agriculture on a higher level [57,59].

#### *5.3. Cultural Heritage under Erosion Threat*

Since the formation of the reservoir in the middle of the 1950s, the confluence of the Kama and Volga Rivers and the left bank tributaries was flooded. As a result, many lower terraces, that were hosting archaeological sites of different periods [60], were completely flooded. The main typology of the sites is presented in Table 2; therefore, out of the total of 1289 sites, 1091 are underwater or totally impacted following the building of the Kuibyshev reservoir. According to their chronology, shown in Table 3, the only Palaeolithic/Mesolithic site that still exists, but is under high threat from coastal erosion, will be further analysed, based on the old Soviet Maps and modern surveys.

Based on the working scenarios regarding the water level increasing and decreasing to 0.5 m and 1 m, respectively it has been observed that increasing the water level, whether, with 0.5 or 1 m, a number of two extra sites will be affected (out of 1091 already underwater or impacted). If we decrease the water level by 0.5 m or1mrespectively, the same number of sites will remain affected—1091. Having such a large surface, water level oscillations do not affect the cultural heritage sites, unless there are variations greater than ±1 m.

#### *5.4. Beganchik Site*

In order to analyse the coastal dynamics of Beganchik site, all the surveys were overlapped, and the site was divided into three sectors (Figure 7), which will be further analysed separately. Beganchik site is located at the mouth of Aktai River, on the second terrace (the first terrace being flooded by Kuibyshev reservoir) of the floodplain which formed before the Holocene [61]; the altitude is between 54–60 m a.s.l. According to the general view of the site (Figure 8a), the northern part of the site is represented by a very steep cliff (Figure 8b) which is continually eroding. Previous preliminary studies [44] have revealed that the erosion rate is about 2–3 m/year.


*Water* **2019** , *11*, 591

**Not affected**

76

 6

 52

 30

 5

 20

 1

 8

 198

**(under**

490

 39

 223

 118

 6

 144

 19

 52

 1091

**water)**

**Figure 7.** Shoreline limit resulted from cartographic analysis and field surveys and the division of the three analysed sectors.

**Figure 8.** (**a**) General view of Beganchik site (drone flight) from 2017; (**b**) Detail over the northern part of the site, fresh parts from the coast are visible in the water (August 2017); (**c**) The change of water colour, due to the clay content of the soil; (**d**) The northern part of the site, where the height of the coast is decreasing.

#### 5.4.1. Sector 1

Sector 1 was not actively eroded between 1958 and 1980 because it was protected by another island (60–90 m north-west, Figure 7), as indicated by the relatively low values of the shoreline retreat (Table 4); in this way, the site was protected from the mechanical action of waves (Figure 8c). Later on, it can be seen that after the island disappeared, the yearly erosion has considerably increased, along with the specific land loss and volume. The direction of the Kama River flow is from north, north-east; being located at the "shelter" of Sector 2 from the speed and currents of the Kama River, this section was in some way protected. However, this sector became likely to be eroded, due to the high erosion rates of Sector 2 and having an elongated shape; this is highlighted of the specific land loss for 1958–1980.


**Table 4.** Detailed morphometric indicators from different observation periods for Sector 1.

Note: \* defines multiplication.

Very high values of the specific land loss, in comparison with other sectors, is due to the height of the coast; which, in some parts can reach 5 m in height. Following the analysis, Sector 1 can be characterised as an extremely dangerous one.

#### 5.4.2. Sector 2

Unlike Sector 1, Sector 2 was and still is under the direct exposure of the Kama River flow and currents. As can be seen in Table 5, the specific land loss is at extremely high rates (Figure 8d). This sector is the most exposed and threatened by erosion. The shoreline retreat is generally stable, varying within 2 m. According to the specific land loss indicator, it can be observed that the destruction occurred, especially within the first two periods, which is typical for the initial stage of lowland reservoir development. During this period, the extremities of this sector are cut off, after which the erosion process stabilises. Between 1958–2008, approximately 70% of the eastern part of the site was eroded. Following that, part of the river's current's strength was redistributed along the north-western part, which explains the sudden decrease in land loss. The height of the coast does not exceed 2 m, therefore, Sector 2 can be classified as moderately dangerous.

**Table 5.** Detailed morphometric indicators from different observation periods for Sector 2.


Note: \* defines multiplication.

#### 5.4.3. Sector 3

Sector 3 is located in the close proximity to Aktai River mouth, where is protected from the mechanical action of waves and Kama River strong currents. This portion of the Beganchik site shoreline is the most stable. As it can be seen from Table 6, there have been no significant changes regarding this part of the coast from 1958 to 2018; except the period 1980–2014, when the specific land loss is higher when compared to other periods, but considerably lower when compared with the other two sectors. The most intensive processes of coastal transformation in the study area were observed in Sectors 1 and 2, open to the destructive effect of the currents and the mechanic action of waves. The erosion intensity may vary from year to year, depending on the water level oscillations in the reservoir. In order to have a more detailed situation on the Beganchik site erosion rates, continuous annual observations are needed.


**Table 6.** Detailed morphometric indicators from different observation periods for Sector 3.

Note: \* defines multiplication.

Particular attention should be paid to Sector 1, in which the most important part of the site is located. If the erosion rates remain stable, the site will be completely impacted in about two or three decades. This imposes urgent mitigation measures from local authorities, along with the sustainable management of cultural heritage sites.

#### **6. Discussions**

Reservoir construction had a significant impact on the flow regime because the current velocity decreased. The currents are very complex, as river flows are under the direct effect of convective flows and wind effects formed in the reservoirs. These are characteristic for Kuibyshev reservoir, where wind effect and bottom relief have a high influence on hydrological conditions [51]. Volga River frames itself into the future increase of global river flow as a consequence of climate change; predictions have shown an increase of 4–8% during 2071–2100 [62]. To be more specific, future trends in the area show an increase in precipitation, temperature and in the use and levels of waters in rivers [63]. The cyclic oscillations have occurred in the Volga River basin in the last half-century; this has influenced the water level in the reservoir, and therefore the erosion rates of the shoreline. The numbers in the Tables 4 and 5 are related to the cyclic oscillations that brought two high-water periods (1951–1962, 1977–1995) and two low-water periods (1963–1976, 1996–present) [51]. For Sectors 1 and 2 high-water levels are associated with low erosion, while the low-water level is associated with higher erosion rates. The research presented in this paper continues our endeavour to monitor the endangered cultural heritage sites from the shoreline of the Kuibyshev reservoir [25–27,60]. Combining old maps with new data collected from field surveys shows high efficiency in establishing the erosion rates of archaeological sites located on shorelines of big reservoirs. When comparing erosion rates with the previous study [25], the average shoreline retreat is close (~3–4 m/year).

#### **7. Conclusions**

In this study, the main changes along the largest reservoir in Europe—Kuibyshev (Russian Federation) were analysed, in strong connection to cultural heritage sites. Following the analysis, Sector 2 has been identified as one with the highest values of width oscillation, from 0.8–1 km to 13–36.8 km. Cultural heritage sites located in the close proximity to big rivers and/or big reservoirs are especially subjected to erosion from water, water level oscillations, and the mechanical action of waves. A diachronic analysis of the archaeological sites located along the Volga River and its main tributaries has highlighted the fact that the most inhabited area was located at the junction of Kama River into the Volga. As highlighted in our analysis, 85% of the cultural heritage around Kuibyshev

reservoir is impacted. However, a more thorough process of monitoring and evaluating the present state of cultural heritage is needed. This has to be done with the cooperation of local authorities and stakeholders. The survey of the only left Palaeolithic site—Beganchik, has shown a fast degradation with no mitigation measures from the local authorities. Beganchik site remains promising for regular rescue archaeological excavations, despite the loss of more than a half of its surface, due to coastal erosion over the last 30 years. Working on scenarios regarding the management of archaeological sites around Kuibyshev reservoir represents one of our future goals.

**Author Contributions:** Conceptualisation, I.C.N. and B.U.; methodology, I.C.N., B.U. and I.G.; software, I.C.N., B.U. and I.G.; formal analysis, I.C.N. and B.U.; investigation, I.C.N., B.U. and I.G.; resources, I.C.N., B.U., I.G. and M.G.; data curation, I.C.N., B.U., I.G. and M.G.; writing—original draft preparation, I.C.N. and B.U.; writing—review and editing, I.C.N., B.U., I.G. and M.G.; visualisation, I.C.N. and B.U.; supervision, I.C.N. and B.U.; review corrections, I.C.N. and B.U.

**Funding:** This work is performed according to the Russian Government Program of Competitive Growth of Kazan Federal University.

**Acknowledgments:** The archaeological database was kindly provided by the Institute of Archaeology of Tatarstan Academy of Sciences (through Leonid Vyazov). James S. Williamson (Memorial University of Newfoundland, Canada) is kindly acknowledged for the English language editing of the manuscript. The authors are grateful for the constructive comments of two anonymous reviewers.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Mathematical Reconstruction of Eroded Beach Ridges at the Ombrone River Delta**

**Irene Mammì 1,2,\*, Lorenzo Rossi 1,2 and Enzo Pranzini 1,3**


Received: 10 September 2019; Accepted: 28 October 2019; Published: 31 October 2019

**Abstract:** Several remotely sensed images, acquired by different sensors on satellite, airplane, and drone, were used to trace the beach ridges pattern present on the delta of the River Ombrone. A more detailed map of these morphologies, than those present in the literature, was obtained, especially at the delta apex, where beach ridges elevation in minor. Beach ridges crests, highlighted through image enhancement using ENVI 4.5 and a DTM based on LiDAR data, were then processed with ArcGIS 9.3 software. Starting from this map, a method to reconstruct beach ridges segments deleted by the transformations of the territory is proposed in this paper. The best crest-lines fitting functions were calculated through interpolation of their points with Curve Expert software, and further extrapolated to reconstruct the ridges morphology where human activity, riverbed migration, or coastal erosion eliminated them. This allowed to reconstruct the ridges pattern also offshore the present delta apex, where the shoreline retreated approximately 900 m in the last 150 years. Results can be further used to implement conceptual and numerical models of delta evolution.

**Keywords:** beach ridges; mathematical reconstruction; curve fitting; Ombrone River Delta

#### **1. Introduction**

Deltas are dynamic depositional landforms, sensitive to changes in both the terrestrial and marine environment [1,2], constantly reshaped by river input, tides, and waves [3,4]. Wave-dominated deltas are characterized by a cuspate morphology, a steep offshore profile, and the presence of more or less curved beach ridges parallel to past shoreline positions [5]. This allowed to reconstruct Holocene delta formation phases for many rivers, and to date, the various morphologies when archaeological remnants were present or if radiometric dating were executed [6,7].

Fluctuations in fluvial sediment supply were also assessed for several deltas, being the most important factor controlling the spatial and temporal variation in beach ridges development [8,9]. Spits are frequently present near the river mouth, sometimes overpassed by the coastal progradation and detectable in the delta area, others eroded by waves. Beach ridges frequently evolve in coastal dunes which reflect a huge amount of sediments that are redistributed by waves and moved inland by winds forming elongated and parallel morphologies behind the coastline [10]. A fast accreting coast makes the formation of several low crests, whereas a slowly prograding or a stable one allows the growth of few high dunes [11,12]. Beach ridges morphology and related pattern could be used as markers of morphodynamic variations: they can give information on past wave regime, climate conditions, sediment supply, and sea level change [13–16].

In historical times, fast delta progradation is frequently induced by fluvial input increase, ascribed to the population development on the catchment area accompanied by intensive deforestation [17,18].

Just short (months) and medium time (years) changes in sediment input can induce delta apex reshaping, resulting in different delta lobes morphology and beach ridges patterns. Fast progradation periods, due to higher sediment input by the river, sees a sharp cusp and the rotation of the river mouth towards the dominant wave front direction, as described for the Ombrone and Arno River deltas [19]; on the contrary, a reduced fluvial input determines a flattening of the cusp up to a rectification of the coastline and downdrift shift of the river mouth [20,21].

Conceptual models of cuspate delta evolution have been proposed by ref. [19,22], whereas other authors ref. [2,21,23,24] used numerical models to explain delta asymmetry and wave reshaping of the river mouth during erosion and accretion phases. For deltas characterized by alternating phases of growth and erosion, not all the beach ridges at the delta apex can be identified, since their most offshore part has been beheaded during the erosion phases [12]; in addition, some segments can be leveled by human activity, especially agriculture, or cut by the river course migration.

A full reconstruction of the delta evolution needs to know the delta shape at its extreme expansions, even if later cut by erosion; this is of outmost importance when conceptual or numerical models based on wings surface are to be developed.

The possibility to observe the secular trends of beach ridges on the Ombrone delta plain represents a valuable opportunity to obtain useful information on its complex geomorphological evolution, where one of the main drawbacks is represented by the loss of geographical information due to natural and anthropic transformations of the territory over time and the lack, or inaccuracy, of ancient cartography. The problem was solved through a mathematical reconstruction of the missing portions of the beach ridges inspired form by the observation of their geometry and how their planar forms remind of specific mathematical functions.

On Earth, a number of geomorphic elements have a regular shape, including small-scale landforms such as cirques, drumlins, sand dunes, alluvial fans, dolines, polygonal soils and cracks, columnar jointing, and scoria cones, to large-scale landforms such as composite and shield volcanoes [25]. Geomorphometry is the science of quantitative land-surface analysis and a variety of landforms are fitted by mathematical functions, often by means of regression analysis [26].

Just to remain in coastal environment, zeta bays in equilibrium with approaching refracted waves are fitted by a logarithmic spiral as proposed by Silvester and Hsu [27], or the equilibrium nearshore profile is mathematically described by the Dean equation [28].

Following a similar approach, a method to reconstruct beach ridges, where they have been eroded or razed by human activity, is hereafter proposed using fitting equations.

Results can be further used for a detailed model on the delta evolution reconstruction, which is not the aim of this study.

#### **2. Study Area**

The Ombrone River delta is in Italy, located on the Southern Tuscany coast (Italy), between Castiglione della Pescaia and Monti dell'Uccellina, in the Maremma Regional Park (Figure 1). The dominant waves direction is from Souh-West [24], and for the present research, is assumed to be a constant wave climate over the historical times [29,30]. The regional longshore transport direction is from south to north [31,32]; subsidence in the alluvial plain was measured in 3 mm/year [33] and this is the only significant vertical movement in the considered time lapse. Sea level rise is about 1 m since the roman period [34].

The ancient gulf of Grosseto started to fill with the material transported by the Ombrone River during the Holocene transgression, but the most internal beach ridge is very likely associated to a sand bar closing a lagoon. When the lagoon was almost silted, the river could directly empty into the sea and started the delta formation. This is dated to the early Roman period [35] and explained with the intense deforestation carried out within the watershed by the Etruscan, first, and by the Romans, later.

During the last twenty centuries, the delta apex grew for approximately 6.5 km at the apex, with a very fast expansion in the 14th–18th century [36], when an intense deforestation occurred. Some authors also see in this delta expansion, the effect of the Little Ice Age [35,37]. Two main erosional phases have been identified, the first after the fall of the Roman Empire, and the second with the Black Death epidemic; in both cases a strong population reduction induced the abandonment of many cultivated areas and the growth of the forest [36].

Furthermore, at the end of the 19th century, these tendencies drastically reversed with the erosion of the coastline at the river mouth, caused by agriculture abandonment, riverbed quarrying, dams and weirs construction; whereas wetland reclamations, through river course diversion, started some centuries before depriving the river load mostly in its finer fraction. Results of these works were poor since sediment compaction was restoring the original wetlands [38].

Since the end of the 19th century, the delta apex started to retreat, and it is now approximately 900 m back from its most prominent position (depicted in a 1881 topographic map). Before shore protection works were carried out in 2013–2014 on the southern lobe (submerged groins) and milder in 2016 on the northern one (short wood permeable groins), the erosion rate was approximately 10 m/year [17].

**Figure 1.** The Ombrone River delta and the position of the two topographic profiles drawn in Figure 3 (Google Earth image).

Delta morphology is slightly different in the two lobes [6,36,37]. The northern side hosts inter-dune swales and, near the river course, some permanent ponds, possibly old channels to connect salterns to the sea. Dunes progressively merge to the north in few parallel ridges; the inner one is attributed to the Etruscan period, the next five to the Roman time, followed by seven Medieval and the last four belonging to the Modern Age [39].

Agriculture activity and urban development at Marina di Grosseto deleted some segments of these dunes. The southern side is almost uninhabited, dunes are closer, and their continuity is more preserved, but at the apex, where they are far lower, seasonal ponds form and in the wetter period most of the apex is flooded. The Monti dell'Uccellina, with their seaward projection, limit sediment dispersion and, even if in the past centuries the prevalent longshore transport was to the north, the delta could symmetrically grow. (Smaller volume but shorter base in the southern side). The River Ombrone delta, hosts many beheaded ridges; but the well-preserved morphology of those remained allows to test a methodology to mathematically reconstruct the lost segments.

#### **3. Material and Methods**

For this study, a detailed beach ridge map was obtained through the combination of several data sets, produced by active and passive remote sensors, and processed with various methodologies. A LiDAR (light detection and ranging) survey acquired in 2008 by the Ministry of the Environment, with 1 × 1 m ground resolution, was used to derive a 3D model of the beach ridges area. Where agriculture flattened the ridges, ad hoc topographic surveys were carried out, and a higher vertical resolution DTM (digital terrain model) was produced from low altitude drone imagery using the structure from motion (SfM) technique.

Satellite images used span from Landsat TM and Landsat ETM+ (from 1986 to 2002); Quickbird (2004); Landsat 8 (2015) and Sentinel 2a (2017); in addition, some Google Earth images were used. Pixel size goes from 2.44 m (Quickbird) to 30 m (Landsat TM).

Remotely sensed images helped in tracing ridge crests where ridges were flattened, and soil moistures helps to discriminate sand from silty/clayey soils present in the interdune swales. This was obtained with specific image processing, such as normalized difference of vegetation index (NDVI), tasseled cup transformation (brightness, greenness, and wetness) [40], principal components analysis (PCA), and RGB (red-green-blue) color composites of the results of the different elaborations. Recent topographic surveys, regional cartography, and old aerial photos acquired in 1944 by RAF (Royal Air Force ) were also used. Images were then wrapped on the Lidar 3D model to assist their photo interpretation.

All these data were processed using ENVI 4.5 (Harris Geospatial Solutions Inc., Broomfield, CO, USA) and ArcGIS 9.3 (ESRI) softwares in order to produce a detailed beach dunes/beach ridges map of the delta area in the WGS84 UTM32 reference system (Figure 2) and leading information was derived from LiDAR data, whose resolution was maintained also after the following image overlapping.

**Figure 2.** Ombrone River delta, beach ridges pattern mapped through the combination of different data sets and elaboration methodologies.

From LiDAR-DTM, several cross-shore profiles were extracted, and minor beach ridges sequence and evolution phases were recognized and relatively dated, starting from the main delta evolution periods described in the literature [36,37].

From the position, geometry, and height of some ridges it was possible to pair them across the gap formed by the river or by agricultural flattering. In Figure 3, profiles crossing the dunes/ridges system on the two delta wings are drawn and ridges associated and numbered (see Figure 1 for profiles location). In the northern profile, all the morphologies are compressed, both because the profile is located far from the delta apex, and because sediment is distributed along a wider area (the northern limit of the cell—Punta delle Rocchette—is 24 km from the river mouth, whereas the southern, one—Monti dell'Uccellina—is 7 km only).

Profiles show the differences in height and distance between ridges, and how some of the highest/largest ones on the delta wings and associated to apex erosion or stability periods, comprise more than one crest, due to an overlapping of more dunes.

**Figure 3.** North (**A**) and south (**B**) profiles from the early 15th (1) century until the limit of the progradation in middle 19th century (9). Due to the difference in deposition rates, horizontal scale cannot be considered a time scale.

Sixteen beach ridges from the early 15th to the middle-19th century have been selected representing different periods of the delta evolution, both when accreting and when eroding. In the first case, beach ridges plant geometry is more curved at the apex with ridges sharply converging (Figure 4); on the contrary, soon after a period of delta erosion their planform is more rectilinear since beach retreats at the apex but progradation goes on at the wings (actually it is a form of beach rotation) [12].

Along each beach ridge, several x,y points were digitized where each ridge crest was more evident; one point every 50 m on average. Along this distance, ridges are almost rectilinear and more points do not increase the accuracy. They were further rotated and shifted in a Cartesian coordinate system with the y-axis coincident with the direction of the terminal river course.

In order to mathematically reconstruct the missing parts of the beach ridges, 33 families of linear and not linear regression models were tested, both for progradation and for erosion beach ridge shapes. The calculation was done using a specific curve fitting software (CurveExpert Pro 1.4, Hyams Development, Chattanooga, TN, USA). The software provides a ranking of the best fitting functions and they are also displayed in the graphing window with statistical results, such as standard error and correlation coefficient. Values x,y for points of interpolated and extrapolated curves are also available and useful for the reconstruction of the missing parts of the delated beach ridges.

#### **4. Results**

Finding equations that well fit the ridges allows to reconstruct their original plan form in the stretches where they have been erased by human activities or by coastal erosion.

Among all the 33 regression models tested, two distinct families of functions, respectively Weibull (Equation (1)) and Morgan–Morgan–Finney (MMF) (Equation (2)), were found out with lowest values of standard error and highest correlation coefficient. The former function better approximates more rectilinear beach ridges of older coastlines or of erosional phases (more flat lines); whereas the latter function better fits curved beach ridges of fast accreting stages (Table 1).

To a lesser extent, Hoerl (Equation (3)) and 3th degree polynomial (Equation (4)) were the following two better equations in the fitting ranking.

$$\text{Weibull}: \text{ y = a - be}^{-\text{cc}^d}. \tag{1}$$

$$\text{MMF} \colon \text{y} = \frac{\text{ab} + \text{cx}^{\text{d}}}{\text{b} + \text{x}^{\text{d}}}. \tag{2}$$

$$\text{Hoeerl} \colon \mathbf{y} = \mathbf{a} \, b^{\mathbf{x}} \mathbf{x}^{\mathbf{c}}.\tag{3}$$

$$\text{Polynomial}: \text{ 3th}-\text{degree } \text{y} = \text{ a} + \text{bx} + \text{cx}^2 + \text{dx}^3. \tag{4}$$

Obviously, the robustness of the method depends on the extension of interpolated beach ridges too, and is correlated to the segment form: linear or arcuate. For very short segments, results exponentially lose significance.

**Table 1.** Relation between beach ridge and best four fitting mathematical functions. ER = apex erosion, PR = apex progradation. In bold, the best standard error for each ridge.


In the table are reported precise year, just to have a temporal reference to be attributed to every mapped beach ridge, but there may be indeed an uncertainty of several years. The correlation coefficients range are from 0.999 to 0.997.

In Figures 4 and 5, the real 1625 (apex erosion phase) and 1750 (apex accretion phase) beach ridges on the south lobe of the delta are drawn together with the best 4 fitting models produced through extrapolation, near the river, of the first 4 km of each ridge. The graphs show Weibull (Equation (1)) better approximates more rectilinear erosional beach ridges and MMF (Equation (2)) curved beach ridges of progradation phases.

**Figure 4.** Real beach ridge 1625 ca., correlated to an erosional period, with the four best fitting mathematical functions.

**Figure 5.** Real beach ridge 1750 ca., correlated to a progradation period, with the four best fitting mathematical functions.

Starting from all the beach ridges mapped in Figure 2, the main evolution phases of the Ombrone River delta, for the considered period, are reported in Figure 6. The considered period is from the early 15th to the middle 19th century, after which reliable topographic surveys are available.

Beach ridges deleted near the river because of watercourse shifting or truncated by a subsequent erosive event have been mathematically reconstructed according to the methodology proposed.

Beach ridges beheaded at the river mouth are represented in cyan and those related to the cuspate progradation in red. In blue and purple, their reconstructed morphologies with mathematical functions.

**Figure 6.** Beach ridges reconstruction for main historical phases from the 15th to the 19th century.

#### **5. Discussion**

Quantitative approach in Geomorphology became established in the mid-20th century [41], but equations fitting coastal landforms arrived a few decades later, like those aforementioned, by Dean [28] and by Silvester and Hsu [27]; both, although developed on actual morphologies, are used to forecast the shape a landform (nearshore profile and zeta-bays) will reach at equilibrium. What is done in this paper is the reconstruction of ancient landforms in those parts which have been later eroded.

In a period in which approximately 31% of the world beaches are eroding [42], and considering that erosion generally starts at the river mouth and gradually expand laterally [22], many present cuspate river deltas lost their apex. Further, agricultural works and urbanization frequently flattened, or at least hid, original morphologies. This makes it impossible to fully reconstruct the natural delta morphology and study its historical evolution.

The proposed method needs high-resolution data, which was obtained from active and passive remote sensing systems. In this way, a detailed reconstruction of beach ridges morphologies was done, and several accretion and erosion phases mapped. Thirty-three different fitting curves were generated though interpolation of control points on the beach ridges/dune crests and their accuracy in describing those morphologies was assessed.

The further extrapolation of the last seaward points allowed to reconstruct the ridge segments which were eroded during the last 150 years, thus obtaining information on the past delta apex shape. On the other sides, interpolation made it possible to identify reaches flattened by agriculture or by other anthropogenic activities. For ridges formed during delta expansions phases (more curved at the tip), the MMF model proved to be more accurate than all the other tested; whereas ridges formed during delta tip erosion, more flattens, were better reconstructed by the Weibull equation.

#### **6. Conclusions**

The mathematical interpolation and extrapolation of the functions fitting the ridge crests at the Ombrone River delta allowed to reconstruct their entire morphologies, not only in the part eroded during the last 150 years, but also where similar processes worked in the Middle and Modern Age.

On these bases, it will be possible to develop more accurate models of delta evolution in order to better explain the different morphology characterizing the updrift and the downdrift delta sides, and to forecast future shoreline displacement in this coast rich in natural and economical elements.

The hoped application of this method to other deltas, using different equation coefficients, could help to solve similar problems and improve the method itself.

**Author Contributions:** Conceptualization, I.M., L.R., E.P.; methodology, I.M., L.R., E.P.; validation, L.R.; formal analysis, I.M.; investigation, I.M.; resources, E.P.; data curation, L.R.; writing—original draft preparation, I.M. and L.R.; writing—review and editing, L.R. and E.P.; supervision, E.P.

**Funding:** This work was supported by the Project Interreg Marittimo-IT-FR Maritime "Management des Risques de l'Erosion cotière et actions de Gouvernance Transfrontalière" (MAREGOT).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
