Key Parameters for Landscape Evolution and Anthropogenisation Estimation in the Kamchia River Downstream Region (Eastern Bulgaria)

Pollen productivity еstimate (PPE) and relevant source area of pollen (RSAP) are critical parameters for quantitative interpretations of pollen data in palaeolandscape and palaeoecological reconstructions, and for analyses of the landscapes evolution and anthropogenisation as well. In light of this, the present paper endeavours to calculate PPE of key plant taxa and to define the RSAP in the Kamchia River Downstream Region (Eastern Bulgaria) in order to use them in landscape simulations and estimations. For the purposes of this research, a dataset of pollen counts from 10 modern pollen samples together with corresponding vegetation data, measured around each sample point in concentric rings, were collected in 2020. Three submodels of the Extended R-Value (ERV) model were used to relate pollen percentages to vegetation composition. Therewith, in order to create a calibrated model, the plant abundance of each pollen type was weighed by distance in GIS environment. The findings led to the conclusion that most of the tree taxa have PPE higher than 1 (ERV3 submodel). Cichoriceae, Fabaceae and Asteraceae have lower PPE.


Introduction
Landscape evolution and anthropogenisation models were numerical models that simulates landscapes changes and human impact over the course of time (Coulthard, 2001). Spore-pollen analysis is one the basic methods used in palaeolandscape reconstructions and estimation of landscape evolution and anthropogenisation (Overpeck et al., 2013). Its accuracy depends largely on the creation of modern calibrated models that reflect the relationship between modern pollen loads and modern vegetation (Sugita, 1994;Vergiev, 2014). Unlike classical qualitative interpretations of pollen spectra, mathematical modeling offers opportunities for relatively accurate quantitative reconstructions of palaeovegetation based on modern dependencies, deriving a number of models based on the pollen-vegetation ratio in modern conditions . This relationship is revealed in the analysis of pollen spectra from modern surface samples located in the investigated area, and the data are compared with the plant communities that produced the corresponding amount of pollen (Andersen, 1970). It is of great importance for obtaining accurate quantitative reconstructions of palaeovegetation in mathematical modeling in GIS environment, to obtain accurate data on key parameters: Pollen Productivity Estimate (PPE) and Relevant Source Area of Pollen (RSAP). Unlike the other parameters required for landscape reconstructions and simulations, these two vary depending on the latitude and the types of plant communities that produced the pollen. This requires that these indicators have to be calculated for each specific landscape and can be subsequently used in modeling (Vergiev, 2014).
Theoretical formulations and research show that the larger sized area collect more pollen (Sugita, 1994). The RSAP reflects the area around the sampling point from which it is possible to produce pollen from the relevant taxa and to assume it to be local. The other part is result of long-distance transmission and its quantity is constant (Sugita, 1994;Vergiev et al., 2014). The aim of the present study was to calculate PPE of 30 key plant taxa and to define the RSAP in the Kamchia River Downstream Region (Eastern Bulgaria). Data will be used to track annual variations in these parameters and will allow to prepare a calibrated model required for palaeolandscape and palaeoecological reconstructions, as well as analysis of the landscapes evolution and anthropogenisation.

2
Materials and methods.

Study area
The Kamchia River downstream region (Eastern Bulgaria) is characterized by diverse semimountainous landscapes and is located through numerous river and sea terraces. It overflows along wide valley and ends with the sand dunes along the longest continuous sandy beach over the Bulgarian Black Sea Coast (Popov & Mishev, 1974). Sea water body proximity determines moderate continental and milder climate characterised with average air temperature of 1.2° C in January and 23° C in July and annual rainfall of 498 mm (Velev, 2002). The vegetation in the investigated area is rich and diverse and characterised by broadleaf deciduous forests and open vegetation communities formed by typical psammophytic and halophytic vegetation (Bondev, 2002).

Spore-pollen analysis
The identification of the pollen-vegetation relationship was carried out on the basis of 10 modern surface samples: 5 surface moss samples (MS) and 5 surface soil samples (SS) (Fig. 1). For analysis, 5 samples were randomly collected from the soil substrate at a depth of up to 5 cm within a sample site measuring 1 × 1 m. The samples were mixed and a quantity of 3 cm 3 was taken from it. Several sub-samples of different types of moss were collected within a limited area of 1 × 1 m and mixed together in one surface moss sample. The laboratory processing of the samples was performed in the Biological Laboratory of the Department of Ecology and Environmental Protection in Technical University of Varna in accordance with the standard acetolysis method of Faegri & Iversen (1989) and Birks & Birks (1980). To determine the pollen spectra, non-permanent glycerol microscopic slides were prepared, in which all pollen grains and spores found in the samples (minimum pollen sum of 250 pollen grains) were identified to the lowest taxonomic level and referred to a certain pollen type.

Phytocoenological investigation
The phytocenological survey of the vegetation around the samples was performed in 2020 following the methodology of Bunting et al. (2013), modified for large areas (Vergiev, 2014) (Fig. 2A). The quantitative participation of the species is described in 4 concentric rings, 8 equidistant transects, starting from the sampling point and directed outwards and evaluated visually in test sites of 1 × 1 m for grass communities and 5 × 5 m for forest and bush communities. The degree of abundance was calculated following the method of Braun-Blanquet scale (1964). The study sites were located according to the methodology of Broström et al. (2004) (Fig. 2B).

Simulations
The obtained distribution and coverage data were digitized in vector format using the QGIS 3.0 Girona software product. The pollen data were processed with the program PolERV v.4.0 (HUMPOL v. 3.1) (Bunting & Middleton, 2005). Three submodels of the Extended R-Value (ERV) model was used to relate pollen percentages to vegetation composition (Parsons and Prentice, 1981;Prentice and Parsons, 1983;Sugita et al., 1999). Pollen productivities and fall speed of different taxa was taken into account (Broström et al., 2004).

Results and discussion
The obtained data on the distribution of modern vegetation are averaged for each studied perimeter and are digitized in vector format in GIS environment. The percentage coverage of each plant taxon was calculated and maps were prepared for each of the samples. For the needs of landscape simulation and testing of mathematical modeling, the studied areas are divided into two main groups: areas that do not produce pollen (water surfaces and urban, urban and bare areas) and areas that produce pollen (mixed oak and hornbeam forests, pastures, meadows, arable land, psammophytic communities and wetlands).
During the spore-pollen analysis of the modern surface samples, the number of pollen grains per 30 key taxa was calculated, and their percentage participation in each pollen spectrum was recalculated to Page | 89 100% (Fig. 3). The limit of 30 taxa is imposed by the pollen sedimentation rate parameter. It is determined experimentally, and the available literature data are for 30 taxa (Broström et al., 2004;Sugita et al., 1999;Fredh et al., 2012).

Fig. 3. Pollen spectra of modern surface samples
The multivariate analysis of the vegetation data (within the different zones of the sampling point) and the percentage participation of each pollen type in the pollen spectra determined the correspondence between the two data sets. Redundancy analysis (RDA) showed a high correlation and a strong relationship between pollen values and vegetation (Fig. 4). Mathematically, this similarity is calculated by the angle of general orientation between the two vectors and the meaning of the two axes.
The permutations of the stochastic Monte Carlo statistical method demonstrate a small angle and a correspondingly low value of p (0.01). Poaceae, Carpinus betulus, Fraxinus, Quercus and Betula show the highest degree of compliance. Large deviations are noted for Artemisia, Plantago lanceolata, Aster-type and Cichoriceae.
RSAP is the distance at which the ratio between the representation in the vegetation and the pollen spectra of all taxa has a linear decreasing dependence, which is also determined by the theoretical formulations in the ERV model (Sugita, 1994). After this distance, this dependence remains linear, but the curve falls into the asymptote and oscillates around one value, so it does not obey the weight factor (Jackson and Lyford, 1999). Using the algorithm set in the program PolERV v.4.0 (HUMPOL v. 3.1), the values of the maximum probability function between the pollen spectra and the vegetation produced by them were calculated. After plotting the XY diagram with respect to the distance from the sampling point, a curve was drawn and the radius of the zone estimated at the distance at which the function approached the asymptote was determined. All three ERV submodels were tested. They give close values, but the radius was the lowest value in the value-balanced logarithmic curve.
The RSAP estimate varies between 4250 m and 4780 m, depending on the submodel used. ERV1 shows the highest value of both the maximum probability function (43760) and the radius -4780 m. With ERV2, the lowest and highest values of the function are between 28475 and 37260. The model shows a small difference, but gives a relatively large radius (4620 m). Submodel 3 shows a relatively regular logarithmic reduction of the curve to reach the asymptote. This model demonstrates a better balance of the curve and shows the smallest radius -4250 m (Fig. 5). Based on the theoretical formulations, 4250 m should be chosen for the radius of the RSAP. This is identical to previous studies for the same region (Vergiev, 2014). Pollen productivity was assessed in 30 key taxa (Fig. 6) using all three ERV sub-models. The best productivity estimate is obtained up to the RSAP distance, however, the average productivity and the standard deviation for each taxon up to 5000 m were calculated. In this way, any change that may occur under the influence of external factors outside the area is ignored. Poaceae was chosen as the reference taxon and all PPE values were calculated relative to it, due to the intermediate relative pollen productivity close to 1 (Sugita, 1994). In the calculations of all three submodels Poaceae has PPE 1.0 and standard error 0.
When using ERV1, six of the taxa had a PPE below that of the reference taxa Poaceae (Fraxinus, Ulmus, Cornus mas, Fabaceae, Asteraceae, and Cichoriceae). All others showed higher values, with Betula (8.59) and Plantago lanceolata (12.03) peaking. Quercus, Plantago lanceolata, Rumex, and Chenopodiaceae show a large standard error, indicating large variations in the PPE of these taxa from the RSAP distance (4250 m) to the maximum study distance (5000 m) and beyond.
ERV2 data show large variations and high standard errors. Only two of the taxa studied had a PPE below those of Poaceae (Asteraceae and Cerealia-type). All others show higher values. Peaks were reported in Juniperus (12.1), Chenopodiaceae (10.07) and Plantago lanceolata (10.1). They also have the largest error, 2.45, 2.89, and 2.86, respectively. Large variations and standard errors define this submodel as the least applicable, as evidenced by the calculation of the RSAP. In ERV3, all tree taxa show values higher than 1.0. Only the grassy Cichoriceae (0.82), Fabaceae (0.67) and Asteraceae (0.34) have lower productivity. Peaks were observed in Plantago lanceolata (11), Betula (9.14) and Fagus (6.76). Compared to other submodels, ERV3 shows the lowest standard error values. This in combination with the regular curve of the logarithmic probability shows this model as the most suitable for carrying out the landscape reconstructions for the investigated area.
In most tree taxa, the values are higher than 1.0, which is in consistence with PPE studies in different parts of Europe (Broström et al., 2004). Herbaceous species that are pollinated by insects generally have low pollen productivity (Real, 1983;Mazier et al., 2008). Such taxon showing low values is Fabaceae.

Conclusions
The radius of the Relevant Source Area of Pollen (RSAP) for the Kamchia River downstream region (Eastern Bulgaria) was estimated at 4250 m, which shows a low degree of severity of the sampling distance factor. Submodel ERV3 demonstrates better curve balance and shows the smallest radius.
Based on the performed spore-pollen analysis of modern surface samples and summarized vegetation data, a linear relationship was established for the pollen-vegetation relationship for 30 plant taxa for the Kamchia River downstream region (Eastern Bulgaria). The obtained data showed that compared to the other submodels, ERV3 shows the lowest values of the standard error in Pollen Productivity Estimate (PPE) and demonstrates a better curve balance and the smallest radius. Therefore, this model is most suitable for palaeolandscape reconstructions and simulations in the Kamchia River downstream region.