Meta-Analysis of Geomorphodynamics in the Western Lower Bakırçay Plain (Aegean Region, Turkey)

: The relation between human activities, climate variability, and geomorphodynamics in the Mediterranean region is widely discussed. For the western lower Bakırçay plain in the ancient Pergamon Micro-Region, geoarchaeological studies have shown changes in geomorphodynamics primarily on a site-basis. We reconstruct past geomorphodynamics in the area based on a meta-analysis of 108 14 C-ages obtained from 25 sediment sequences mainly from colluvial and alluvial deposits by analyzing cumulative probability functions of the 14 C-ages. Accounting for biases in the database, we applied different approaches and compared the empirical probability functions with simulated functions. Reconstructed geomorphodynamics in the western lower Bakırçay plain during the Holocene principally coincide with a trend of climate-driven sensitivity to erosion and population dynamics in the eastern Mediterranean, but are also related to the local settlement history. Our data analysis shows that transformations of the Pergamon Micro-Region between the Hellenistic and Roman Imperial times is contemporary to increasing geomorphodynamics that peak in Roman Imperial times. However, a cause–effect relationship between geomorphodynamics and settlement dynamics should be further evaluated. A comparison with data from other settlement centers in Anatolia shows that a coincidence between the peak in geomorphodynamics and a peak in settlement activity are not obvious and may be influenced by soil conservation measures, preferred settlement location, and inherited soil exhaustion. the palaeogeography of the area. Coastal dynamics had an inﬂuence on settlement conditions; a change from progradation to aggradation in the surroundings of the a settlement hill in the delta area is observed.


Introduction
Ever since Claudio Vita-Finzi's seminal book The Mediterranean valleys: geological changes in historical times [1], the varying importance of climate and human control on soil erosion by water in the Mediterranean region has been discussed [2][3][4][5][6]. In a review on Holocene environmental change in the eastern Mediterranean, Dusar et al. [7] identified several phases in which sediment dynamics changed. They stress the decreasing importance of climate control over erosion sensitivity during the Holocene and the coincidentally increasing importance of human impact [7]. The authors additionally state that both erosion sensitivity and sediment dynamics mainly increased from the Early Bronze Age to Late Antiquity. In addition to a general trend of changes in the importance of human impact on sediment dynamics in the eastern Mediterranean, local and regional studies offer a nuanced view on the topic [8]. This becomes especially important when considering different types of archives along a sediment cascade [8][9][10][11][12].

Sediment Dynamics in the Western Lower Bakırçay Plain
For Elaia, the "maritime satellite" of the ancient city of Pergamon [34,35] (the modern city of Bergama) (Figure 1), geoarchaeological and sedimentological studies are also available [36][37][38][39][40][41]. In the bay of Elaia, the construction of breakwaters and moles during antiquity had an impact on local sedimentation rates and styles [39,40]. This change in sediment dynamics was, however, not only due to the constructions of maritime structures and therefore a change in the local depositional environment, but also due to human impact in the contributing catchments [39]. The increasing presence of Glomus-type fungi in the sediments from the harbor around 2.7 ka BP, along with an increasing sedimentation rate, indicators of torrential floods, and palynological proxies, point to increased land clearing and subsequent erosion inland after 2.3 ka BP [37,41]. This impact was much more pronounced in the vicinity of the major settlement around the bay compared to more distal locations [37]. The occurrence of torrential floods and consequent siltation presumably resulted in the abandonment of several coastal settlements and made the harbor inaccessible for certain types of ships whose draught exceeded water depth [38,39]. In addition, the development of the Madra delta plain in the direct vicinity of the western lower Bakırçay plain is studied from a sedimentological and geoarchaeological point of view [42][43][44]. Due to the relatively low chronological resolution, these studies do not give hint to changing geomorphodynamics. However, the studies clearly show the coastal development and the palaeogeography of the area. Coastal dynamics had an influence on settlement conditions; a change from progradation to aggradation in the surroundings of the a settlement hill in the delta area is observed. Figure 1. Overview map of the western lower Bakırçay plain. Sampling locations of the sediment sequences referred to in the current study are highlighted; main grabens are named and selected rivers/creeks and settlements are shown (some are indicated by bold numbers). For details on the sediment sequences see Table 1. Database: TanDEM-X digital elevation model [45,46], rivers partly digitized from Open Street Map-data [47]; locations of sediments according to original publications [37][38][39][40][41][48][49][50][51].
Studies from the western lower Bakırçay plain revealed changed sediment dynamics in the context of historical settlements [37,49,50]: Fan development originating from the slopes around the acropolis of Atarneus started in the Bronze Age, around 3.1 ka BP. Colluvial sediments in the area date to the Middle-Late Iron Age (ca. 2.7 to 2.4 ka BP.) [49]. Despite a clear change in human activities around Atarneus in the last 3 ka, sedimentation processes did not change much. Schneider et al. [49] attribute this to the construction of terraces on the slopes around Atarneus that reduced soil erosion rates even under intensified land use. In the surrounding of the Late Chalcolithic-Early Bronze Age settlement at Yeni Yeldegirmentepe, fluvial sediments are intercalated by sediments showing a clear anthropogenic imprint (archaeological remains were uncovered from the respective layers); a change in sedimentation rate coincides with early activities at the site ( Figure 3 in [50]). In contrast, there is no change in sediment dynamics related to pre-modern human impact recorded for the Geyikli valley and the environs of the archaeological site near Sultantepe [48,51].

Objectives and Outline
In addition to the study from the Geyikli valley that did not recover any traces of human impact on sediment dynamics [48], all cited studies from the western lower Bakırçay plain are directly related to human settlements of a specific period, such as the Graeco-Roman harbor at Elaia [39,41] or the Chalcolithic-Bronze Age settlement at Yeni Yeldegirmentepe [50]. A meta-analysis synthesizing the detailed analyses for the western lower Bakırçay plain is not yet available. Therefore, the main objective of the current study is to re-evaluate the existing sedimentological and geochronological data to reconstruct the general trend in sediment dynamics in the western lower Bakırçay plain.
Our meta-analysis is based on cumulative probability functions (CPFs) of radiocarbon ages sampled from sediment sequences obtained from the western lower Bakırçay plain and its intermediate vicinity. We compare the CPF of the observed ages with a null model of the cumulative probability that can be expected from the ages randomly distributed over the available sediment sequences. The approach is derived from archaeological demography studies [52,53]. Along with an estimation of sedimentation rates and a Bayesian chronological model of facies change, we use the null model approach to reduce the impacts of depositional and sampling biases on the identification of phases of increased geomorphodynamics (see Section 2.3.2).
The generic term geomorphodynamics is used in the current study for two reasons. First, the term includes various surface processes-different processes of erosion, deposition, and reworking are covered. With our approach, we analyze the sediments related to these processes. Since different processes, e.g., increased erosivity or reduced vegetation cover due to human impact, can result in the same sedimentary signal, a term covering various processes is most appropriate. This is especially important as climate and human triggers are entangled [54] and one trigger can level or accelerate another one [4,55]. The triggers of changed geomorphodynamics are assessed by comparing the identified phases of increased geomorphodynamics with different proxies of climatic variability, vegetation development, and population dynamics from Anatolia, southern Greece, and the Balkans. Second, the term geomorphodynamics is used if a change in the sediment sequences is indicated by different proxies, i.e., either combination of the following three: a change in the cumulative probability of 14 C-ages, a change in the sedimentation rate, or a facies change.

Study Area: The Western Lower Bakırçay Plain
The western lower Bakırçay plain is the lowest of the three major plains of the river Bakırçay and lies south-southwest of the Bergama fan ( Figure 1). For practical reasons, we included the coastal areas close to Dikili and Zeytindag to the western lower Bakırçay plain, although this is topographically not consistent. The Bakırçay rises in the Ömer Dagı, flows mainly in east-southeastern direction and drains a catchment area of circa 3350 km 2 into the Aegean Sea [56,57]. The western lower Bakırçay plain covers an area of around 140 km 2 .
The climate in the region is a typical Mediterranean hot and dry summer climate of the temperate regions (following the Köppen-Geiger classification [58,59]). Annual rainfall averages 636 mm (Dikili) or 711 mm (Bergama) with a maximum monthly average in December (115 mm, Dikili) and a minimum monthly average in July (<10 mm). The annual average temperature is around 16°C [60].
The development of the western lower Bakırçay plain is mainly driven by the formation of a horst and graben structure during the Middle Miocene to Early Pliocene (Bergama, Zeytindag, and Altınova-Lesvos grabens and the Kozak and Maruflar horsts) [61]. Bedrock consists of mainly volcanics (pyroclastic and andesitic rocks in the north and east), intrusive rocks (granodiorite in the headwaters of the north), and carbonate and continental clastic rocks (in the south) [62]. During the Pleistocene, pediments and relatively large fans (such as the Bergama fan) developed on the southern margin of the Kozak horst. Pediments are mainly located in the northern part of the western lower Bakırçay plain due to the asymmetric shape of the horst-and-graben structure; average slope and size of the catchments of creeks draining the Kozak mountains to the northern part of the western lower Bakırçay plain are greater compared to the catchments of creeks draining the Yunt Dagı mountains in the southern and eastern part of the plain [60,61]. Several andesitic hills (Turkish tepe, e.g., Sultantepe and Kalerga Tepe) occur on the northern edge of the western lower Bakırçay plain ( Figure 1). A detailed geomorphological characterization of the entire Bakırçay catchment and the Madra catchment, including the adjacent coastal areas, is given in the contribution of Yang et al. in the current special issue [60].

Database
We collected n = 108 14 C-ages from s = 25 sediment sequences obtained from the western lower Bakırçay plain and the surrounding areas ( Figure 1). The sampling locations comprise alluvial terraces of the Geyikli drainage basin in the north-western part of the western lower Bakırçay plain; the foot slopes of the small volcanic Yeni Yeldegirmintepe; the transition zone between the piedmont of the Kozak horst and the plain in the northern part of the western lower Bakırçay plain; and the northern part of the Bay of Elaia, where sediments mainly originate from the slopes and small drainage basins of the Yunt Dagı mountains ( Table 1). The archives at these different sample locations react differently on disturbances such as land-use change or changing precipitation. Deposits from low order catchments resulting in, e.g., colluvial layers, are directly related to local human impacts if there is e.g., a settlement in the close vicinity [10,11]. Sediment archives from high order catchments, e.g., the Bakırçy floodplain, by contrast, react less directly to local human impact-their signal rather gives an average of changes within the wider catchment area. Compared to colluvial deposits, floodplain deposits of larger catchments are thus less sensitive to local human impact of the same scale, but can react equally sensitive to changes in precipitation. In addition, the residence time of the deposits within these archives differs considerably [63].
For all available 14 C-ages (see Figure 2), metadata were collected from the original publications, i.e., sampling depth, total depth of the sequence, uncalibrated 14 C-ages, lab errors, and information on a potential reservoir effect. Additionally, the lithostratigraphy and physical/chemical data of the sediment as given in the original publication were used to categorize the data into change classes [64][65][66][67] ( Figure 2). Ages were classified as "change-before ages" if a major change in sediment facies is indicated in ≤25 cm below the sampling depth (n = 37) and as "change-after ages" if a facies change is indicated in ≤25 cm above the sampling depth (n = 23). Ages were classified as "no-change ages" if they are stratigraphically not related to changes in sediment facies (n = 27). For n = 21 ages, no information on a facies change is available. For all sediment sequences with available 14 C-ages, we also collected lithostratigraphic information from the original publications. As a result that the terms used to describe the lithostratigraphic units differ between the publications, we harmonized the data. The 14 C-ages were on a first level categorized as ages from terrestrial or marine sedimentary units ( Figure 2). Colors indicate different change-ages (black = no change, green = change-before, red = change-after) and ages where no stratigraphic information are available (grey). Horizontal dashed lines separate different sediment sequences. Probabilities of the dates are not normalized to sum to unity. Secondary y-axis: Location codes as given in Table 1 and Figure 1. Ages obtained from [37][38][39][40][41][48][49][50][51]. Table 1. Overview of the sediment sequences included in our meta-analysis. More details can be found in the Supplementary Materials and the original publications. The locations are shown in Figure 1. For n = 21 samples, no detailed information on sediment sequences and geomorphological setting is available. s = number of sediment sequences from location; n = number of 14 C-ages from location.
Maps were created in QGIS Version 3.4.11-Madeira (2019-08-16) [75]; maps and figures were post-processed in Inkscape Version 0.92.4 (14 January 2019) [76]. A detailed description of the computational procedure can be found in the Supplementary Materials.

Calibration and Observed Cumulative Probability Function (CPF)
All 14 C-ages were calibrated in R using the rcarbon package. The used calibration curves are IntCal13 for terrestrial ages and Marine13 for marine ages that are affected by a reservoir effect [77]. Following Crema and Bevan [53] (see [78]), we did not normalize the age probability distributions of the calibrated ages to unity to avoid artificial peaks in the observed CPF ( Figure 3a). The observed cumulative probability function is calculated by summing the dating likelihood of all non-normalized calibrated 14 C-ages for all ages dating between 0 BP and 11,700 BP (so covering the Holocene Epoch). The formal subdivision of the Holocene is based on Walker et al. [79]: Boundary between Pleistocene and Holocene at 11.7 ka BP; between Early and Middle Holocene at 8.2 ka BP; and between Middle and Late Holocene at 4.2 ka BP.
The main principle behind the calculation of CPFs from sediment archives is the assumption that an increase of the CPF is proportional to increased deposition when using data from several sediment sequences (e.g., [80]); the record of a single sequence is believed to be fragmented [81,82].

Biases Affecting Cumulative Probability Functions
The application of cumulative probability functions (CPFs) for sedimentological analysis ignited debate (e.g., [80,83,84]; see also [85]). The main critique of the application of CPFs to reconstruct sediment dynamics and especially fluvial activity arises from the nature of depositional and taphonomic processes and the sampling of the dating material. Several biases hamper a straightforward relationship between CPFs and sediment dynamics.
One issue is the effect of the shape of the calibration curve on a CPF; steeper parts of the curve cause peaks in the CPF [86]. Furthermore, the formation of dating material (especially charcoal and plant remains) is not necessarily contemporaneous to transport or final deposition of the material and, thus, the processes related to sediment dynamics [80]. Additionally, in sedimentary archives, dating material and the sediment layers containing the samples are susceptible to erosion and reworking. Thus, preservation cannot be premised [80,82,85,87,88]. The preservation potential of dating material and the sediment layer containing the dating materials decreases with increasing time difference to the moment of sampling (similar to the Sadler-effect [89]; see, e.g., [90,91]).
A 'researcher' bias may cause further problems. For instance, in geoarchaeological studies focusing on a specific period, researchers may try to well capture this period with their geochronological framework [92,93]. Furthermore, the uncovered depth of a sediment sequence may lead to an additionally increased likelihood that younger samples are selected; the thickness of the sediments may be larger than the sampling depth. The availability of dating material is further affected by a 'production' bias: Human-induced (wild) fires in a specific period may increase the likelihood that dating material from this period is sampled, although sediment dynamics may not have changed [85]. Sample size is crucial when dealing with biases in the record of 14 C-ages [94]. In small data sets, single ages have a greater impact on the shape of the curve of a CPF.
To cope with potential biases in our data set, we estimated a 'null model' of simulated CPFs and compared them with the observed CPF to test for the reliability of the periods of increased sediment dynamics in our data set from the western lower Bakırçay plain. Additionally, we did a CPF-based estimation of sedimentation rates [8].

A Null Model
A common approach in archaeological demography to reduce the effect of biases is to compare the empirical data from 14 C-ages with a null model, i.e., a hypothesis of the shape of the CPF [52,92,[95][96][97]. Our null model is based on the expected likelihood that a point in time is covered with the given sediment sequences. Therefore, we estimated the period covered by each sediment sequence. The underlying hypothesis is that the change in the cumulative probability increases with time, e.g., due to the higher probability that a piece of dating material is preserved, uncovered, and sampled. The minimum age of a sediment sequence s is set to 0 BP; the maximum age is estimated as a function of the total thickness of the sediment sequence and the average sedimentation rate given by the depth and the median age of the calibrated 14 C-ages of the sediment sequences. Sequences for which only one 14 C-age is available, the maximum age of the sequence is estimated by dividing the thickness of the sequence by the sedimentation rate of the given 14 If more than one sample is available for a sequence, the maximum age of the sequence is estimated based on a power regression age-depth model. Although residuals of the model may be high, this approach is appropriate for the estimation of the maximum age and to account for the preservation bias of the sediment record. Other approaches, such as the estimation of the mean sedimentation rate of a sediment sequence [8,86], are less sensitive to changes in sedimentation rate with increasing range between sampling and recorded age (the so-called 'Sadler effect' [89], see [8]).
In total, 1000 CPFs were simulated. For each of the simulated CPFs, several calendar years were randomly sampled from the period between 0 BP and 11,700 BP. The likelihood that a year in that period is sampled is equal to the proportion of all available sediment sequences covering the respective year. The number of samples taken for each simulated CPF is equal to the number of available 14 C-ages from the study area. With this approach, the sampling density in the original data set is taken into account. For the simulation of a CPF, the random calendar ages where back-calibrated to raw 14 C-ages and again calibrated (the estimated error is 40 years for all ages regardless the age of the sample). This approach is necessary to ensure that the observed and simulated CPF both have the same shape for same ages.
The variability in the simulated CPFs is considered by calculating confidence intervals of the simulated CPFs. We furthermore subtracted separately each simulated CPF from the observed CPF. The likelihood that the observed CPF is not increased in a year due to variation in the calibration curve or a sampling bias is equal to the proportion of simulation runs in which the observed CPF is higher than a simulated CPF (the null model). We also normalized the observed CPF by dividing through the cumulative probability of all simulated CPFs. This approach is similar to normalization by a CPF of equally distributed ages [86], but also takes the likelihood of sampling into account.

CPF-Based 'Sedimentation Rate'
To account for further sampling biases (e.g., samples from a specific period preferably sampled), we also took the depth of the 14 C-ages into account. Therefore, we randomly sampled depths from all sediment sequences and simulated a possible age of the sediments a certain depth. The age probability function is a function of the overlying and underlying age in the sediment sequence. We sampled a calendar year from the age distribution of both the over-and underlying sample. Subsequently, one age is drawn from the range between both ages. The process is repeated 1000 times to get an age probability distribution for the random depth. Age distributions were calculated for 1000 random depths. From all simulated ages, a kernel density estimate is computed following the procedure proposed by Bevan and Crema [53]. This kernel density of ages is proportional to sedimentation rates. Additionally, the number of ages sampled in one simulation run is equal to the number of empirical 14 C-ages to calculate confidence intervals of the kernel densities and to take the original sample size into consideration.

Facies Change
At the random depths described above (Section 2.3.4), also the facies information given in the original publication of the sediment sequences was recorded. Thus, for all the different facies units, kernel densities of the possible ages were calculated.

Cumulative Probability Functions
Cumulative probability functions of the available 14 C-ages from the western lower Bakırçay plain are depicted in Figure 3a-d.

Normalization
Except for a period between 2.3 ka BP and 2.8 ka BP, the CPFs of normalized and non-normalized 14 C-ages do not show major differences (Figure 3a). Around 2.7-2.8 ka BP, a distinct peak in the CPF of the normalized data can be observed that is not visible in the CPF of the non-normalized data. In the period after the peak, the cumulative probability of the non-normalized data is in contrast higher than the cumulative probability of the normalized data. This difference is due to a pronounced wiggle in the IntCal13 calibration curve around 2700-2500 14 C-ages BP prior to the so-called Hallstatt-plateau (e.g., [98]). In the following, only non-normalized ages are described.

General Trends in the CPF
The CPF of all ages shows several changes in the cumulative probability, starting at 10.1 ka BP (Figure 3a). For several phases, the dating probability is zero (9.5-9.2 ka BP, 8.2-8.0 ka BP, 5.9-5.3 ka BP, and 1.4-1.3 ka BP). Four major plateaus of increased cumulative probability are recorded (9.0-8.6 ka BP, 5.3-5.0 ka BP, 2.8-2.4 ka BP, and 2.0-1.6 ka BP). After 3.0 ka BP, the cumulative probability increases more or less continuously, reaching a maximum at around 1.7-1.6 ka BP. There is a marked depression in the CPF between 1.6 and 1.3 ka BP that is followed by a more or less continuous record until present.

CPF of Terrestrial and Marine Ages
The CPFs of 14 C-ages sampled from terrestrial sediment layers and the CPF of 14 C-ages sampled from marine ages differ clearly (Figure 3b). This is especially true for the several peaks in the terrestrial CPF before 3 ka BP. While the earliest increase in the CPF of terrestrial ages is recorded for 10.1 ka BP, the CPF of marine ages first increases at 7.8 ka BP. Major phases of an increased CPF of terrestrial ages, viz. 9.0-8.6 ka BP, 5.3-5.0 ka BP, 4.6-4.4 ka BP, and 3.6-3.2 ka BP, are not equally evident in the marine record. In contrast, the CPF of the ages from marine sediment layers is higher than the terrestrial CPF between 7.8 and 7.0 ka BP and between 6.2 and 6.0 ka BP. After 3 ka BP, the increase of the cumulative probability occurs more or less simultaneously in terrestrial and marine records. The CPF of the marine ages decreases after 2.0 ka BP and again peaks at around 1.7 ka BP. During this period, the cumulative probability of the terrestrial ages is continuously high. The cumulative probability of marine ages decreases to zero after 1.4 ka BP. For the following analysis, only the terrestrial ages are considered.

Phases of Increased Terrestrial Cumulative Probability
Statistically, nine phases of increased cumulative probability of terrestrial ages were detected (Figure 3c These phases all lasted more than 100 years; phases with an interruption of less than 100 years are aggregated. In these phases i-ix, the cumulative probability increased the 1000-years running mean of the cumulative probability (and at least the probability of three ages contributed to the cumulative probability during that phase). The running mean shows that the phase between 9.0 and 8.6 ka BP is followed by a super-ordinate phase of low mean cumulative probability, which lasted until 5.3 ka BP. On average, the cumulative probability tends to increase between 5.3 ka BP and 1.6 ka BP, but is interrupted by a phase of low mean cumulative probability between 3.8 ka BP and 2.8 ka BP.

Comparison with the Null Model
Most of the nine phases of increased cumulative probability are statistically significantly different from the null model (simulated CPFs) at 67% confidence (Figure 3d, yellow envelope and light green columns). Three phases are statistically significantly different from the null model at 95% confidence, viz the phase between 9.0 and 8.6 ka BP, the phase between 4.1-3.9 ka BP, and 2.0-1.6 ka BP (Figure 3d, dark green columns). For several time intervals, the cumulative probability is statistically significantly increased, but cumulative probabilities were calculated from the dating probabilities of less than 3 ages or are not higher than the 1000 years running average (e.g., at 8.4-8.3 ka BP or between 7.0 and 6.8 ka BP).

Cumulative Probability Functions of Change Ages
In total, 35 terrestrial change-ages are available from the sediment sequences, indicating change after (n = 15) or change before (n = 20) the given age (Figure 3e). The CPFs of the change-after and change-before ages point to several phases of variation in deposition style; before 4.8 ka BP, the CPF of change-ages is not congruent to the CPF of all terrestrial ages but is especially similar after 2.8 ka BP. No-change ages are recorded before 5.4 ka BP; a single change-age dates between 5.4 and 5.0 ka BP.
A phase of change-after ages is recorded for between 4.8 ka BP and 4.2 ka BP, transitioning into a phase of change-before ages between 4.4 ka BP and 3.9 ka BP. Thus, the change-ages point to a change in sedimentation style between 4.4 ka BP and 4.1 ka BP. A phase of overlapping change-before and change-after between 3.7 ka BP and 3.4 ka BP does not point to a specific phase of change. Between 2.9 ka BP and 2.4 ka BP, only the cumulative probability of change-before ages is increased. The cumulative probability of change-after ages increases after 2.5 ka BP. Between 2.5 and 1.7 ka BP, the cumulative probability of change-after ages is higher than the cumulative probability of change-before ages. Between 1.7 and 1.5 ka BP, only change-before ages are recorded, followed by a phase of no-change ages. The latest marked phase of change-ages is recorded between 1.4 and 1.0 ka BP (change-after ages between 1.4 and 1.1 ka BP and a change-before ages between 1.1 and 1.0 ka BP.)

Sedimentation Rate
The sedimentation rate calculated based on CPFs of estimated ages at random depths points to five distinct phases of higher and lower sedimentation dynamics compared to the null model (Figure 3f). The modeled sedimentation rate is relatively constant at a low level between 10.1 ka BP and 5.4 ka BP. During this period, the mean simulated CPFs of the age-depth model does not exceed the null model (CPF of ages randomly distributed over the sediment sequences) except of the period between 7.5 and 7.1 ka BP, where the 95% confidence interval of the simulated ages increases the confidence interval of the random ages.
The 95% confidence interval of the simulated ages also increases the confidence interval of the null model between 5.1 and 4.3 ka BP and between 2.9 ka BP and present, indicating increased sedimentation rates. Between 4.2 and 1.7 ka BP, the simulated sedimentation rate constantly increased; the average of the simulated CPF exceeds the null model between 2.3 ka BP and 1.7 ka BP. After 1.5 ka BP, the sedimentation rate is not increased.

Facies Change
Several phases of a marked change in the terrestrial facies of the sediment sequences can be identified (Figure 4). The most pronounced change in facies occurred around 7.1 ka BP; being, however, only covered by a small number of dates and possibly due to terminological differences between the original publications. The decreasing proportion of overbank and channel-related sediments and an increase in fluvial sediments might not necessarily point to an absolute change in the depositional environment. The changes in the facies composition at around 6.3 ka BP and 5.8 ka BP are characterized by a decrease of sediments indicating coastline propagation. At around 4.2 ka BP, the proportion of sediments forming alluvial fans and sediments indicating a progradational depositional environment increases. After 2.7 ka BP, the proportion of colluvial sediment layers and overbank sediments increases. At round 2.0 ka BP the summed contribution of overbank, colluvial, and fan sediments reaches a local maximum; the contribution of channel-related sediments increases afterwards. The composition remains relatively stable with minor variations after 1.6 ka BP.

Early Holocene (ca. 11.7-8.2 ka BP)-Aceramic Neolithic
Highest climatic erosion susceptibility in the eastern Mediterranean region after the Pleistocene presumably occurred during Early Holocene peaking around 10 ka BP (Figure 5m) [7], coinciding with an Early Holocene wet period in the Balkans and Anatolia (Figure 5n,o) [99]. During this period, rainfall erosivity was relatively high. This effect may be coupled with a relatively low vegetation cover after the dry late Glacial [100][101][102]. Notwithstanding the high erosion susceptibility during the Early Holocene, sediment dynamics were relatively low and mainly triggered by climate variability and concurrent environmental change. Human impact during Early Holocene was neglectable but increasing [103]. The distribution of farming and a west-ward human migration from the Fertile Crescent associated with the onset of Neolithisation caused an increase in population in Anatolia and the Balkans (Figure 5a,d), however without reaching the Aegean region before ca. 8.5 ka BP [103].  [105]; d: [104]; e: [106]; f, g, h: [105]; i, j: [106]; k: [107]; l: [41], interpolated; m: [7]; n, o: [99]; p: [7]; q: [8]; r-u: this study.
Our meta-analysis from the western lower Bakırçay plain shows one clear, statistically significant increase of the cumulative probability of 14 C-ages during the Early Holocene, ranging from 9.0 to 8.6 ka BP ( Figure 5r) and thus coinciding with the Ceramic Neolithic [108,109]. As a result of lacking metadata for some of the 14 C-ages in our data set, we could not reliably estimate the sediment accumulation in the western lower Bakırçay plain for the Early Holocene. Interpretations have therefore to be taken with caution.
The increased cumulative probability during that period in the western lower Bakırçay plain fit to the general trend of increased population dynamics and vegetation change in southeast Europe and Anatolia. In southern Greece, more clearly than in (Southern) Anatolia, palynological data show an increase in anthropogenic pollen indicator and grazing indicators (Figure 5g-i, [105,106]). Archaeological evidence points to the earliest known human activities around the western lower Bakırçay plain possibly between the 7th and the 5th millennium BCE (Late Ceramic Neolithic to Chalcolithic), although the settlement history might go back to the Aceramic Neolithic [110][111][112]. The prehistoric sites in the western lower Bakırçay plain are limited to the northern part of the plain and the Gümüş valley. The settlement cluster in the Gümüş valley ( Figure 1) developed in the 5th and 4th millennium BCE [112]. The sediment sequences included in our data set are not located in the vicinity of these sites. Assuming a localized impact of the prehistoric sites, the contribution of an anthropogenic trigger to the increased cumulative probability before 8.6 ka BP is still unclear. Potentially above-average wetness (Figure 5n,o) during the phase may have also contributed to increased geomorphodynamics (cumulative probability of 14 C-ages). Berger et al. [113] show an increase in humidity during the period of increased geomorphodynamics that goes along with increased population dynamics. In studies on sediment dynamics in other regions of the eastern Mediterranean [7], the Early Holocene is relatively seldom covered. However, Fuchs et al. [9] reconstructed phases of increased sedimentation rates during the Neolithic on the Peloponnese peninsula (S-Greece); based on this data, Fuchs [55] concludes that Early Holocene sediment dynamics are related to both climate and human impact (i.e., high rainfall reconstructed form δ 18 O values; archaeologically testified increase in settlement dynamics). Studies on the Çarşamba alluvial fan (Konya Basin, Central Anatolia), where the important site of Çatalhöyük is located, revealed an early human impact. Additionally around Sagalassos, a distinct peak in the CPF of alluvial sediments occurs [8], which coincides with the peak in cumulative probability in the western lower Bakırçay plain. Phases of increased soil erosion around Arslantepe (eastern Anatolia) partly overlap with the period reconstructed for the western lower Bakırçay plain. At Arslantepe, the increased geomorphodynamics occurred before the known onset of the settlement period; thus, also here climate might have triggered erosion [114].

Middle Holocene (ca. 8.2 to 4.2 ka BP)-Ceramic Neolithic to Late Chalcolithic
For the eastern Mediterranean region, the Middle Holocene can be characterized as a phase of more or less continuous decline or stability of the climatic erosion sensitivity [7]. In the Balkans, reconstructed hydroclimate followed a trend of continuous aridization from the wettest Holocene period around 8.0 ka BP [99] (Figure 5n) until ca. 4.9 ka BP.
In the area of modern Turkey, the aridization trend is less pronounced, but culminated in the so-called 4.2-ka-BP-drought-event [99] (Figure 5o) that dates to a period around 4.3/4.2-3.9/3.8 ka BP [115,116]. Estimated sediment dynamics in the eastern Mediterranean during the Mid-Holocene remained stable at a low level, but show a phase of increased dynamics between 5.5 and 5.0 ka BP [7] (Figure 5p). Human activities in Anatolia and Greece reconstructed from archaeological cumulative 14 C-ages and the aoristic sum of site counts (Figure 5a-e) first declined during the Mid-Holocene, but increased after 5.2 ka BP (Late Chalcolithic-Early Bronze Age), when also the estimated sedimentation dynamics increased [7] ( Figure 5p).
The cumulative probability of 14 C-ages from sediment sequences and the reconstructed sedimentation in the western lower Bakırçay plain also increased during the aforementioned period (Figures 3 and 5r-u). An earlier peak in sediment dynamics occurred in the (Late) Chalcolithic (especially after 5.4 ka BP), followed by a more pronounced average increase in sediment dynamics during the Early Bronze Age. This increase is in accordance with the model for the entire eastern Mediterranean and also coincides with general changes in vegetation cover in the area of modern Turkey and Greece (Figure 5f-k). In the Aegean region of Anatolia, settlement activities also increased during the Late Chalcolithic; the number of archaeological sites in the area of modern Turkey increased around 5.2 ka BP. The palynological record originating from the region shows an increase of anthropogenic pollen indicators (including a relative increase of OJCV-pollen and grazing indicator species) between 5.7 ka BP and 4.6 ka BP. The palynological record from the Bay of Elaia [41] close to the Bakırçay delta shows a reduction of forest cover during the period of increased geomorphodynamics.
Although an anthropogenic driver of increasing geomorphodynamics in the late Middle Holocene appears to be clear given the fact that human activities increased and palynological data indicates increased land-use pressure, changing hydro-climatic conditions (i.e., increased erosivity) cannot be excluded as a trigger. In the Balkans, a wet phase occurred between 4.6 and circa 4.0 ka BP [99] ( Figure 5n). However, the hydro-climatic reconstructions from the area of modern Turkey indicate a wet phase between 5.4 and 4.4 ka BP [99] (Figure 5o).
In the western lower Bakırçay plain, the number of settlements increased from the Early-Middle Chalcolithic to the Late Chalcolithic-Early Bronze Age [110]. During this period the western lower Bakırçay plain was the first time continuously settled. This also implies that the increased geomorphodynamics were triggered by settlement activities. To what extent the reconstructed geomorphodynamic presented here are representative for the whole western lower Bakırçay plain has to be further surveyed. Some of the sediment sequences in our data set were obtained from locations in the vicinity of Late Chalcolithic-Early Bronze Age settlements, such as Elaia, Yeni Yeldegirmintepe, and Başantepe. Neolithic sites were, however, not found in the vicinity of the locations of the sediment sequences. Therefore, not only the absolute increase in sites during the Late Chalcolithic but also the site pattern may have had an effect on the observed increase in geomorphodynamics.

Late Holocene (4.2 ka BP to Present)-Early Bronze Age to Post-Classic Period
According to the review analysis of Dusar et al. [7], the Late Holocene in the eastern Mediterranean is characterized by low climatic sensitivity to erosion, but high sediment dynamics. Thus, increased soil erosion and a subsequently increased deposition in the lowlands is not due to higher rainfall erosivity, but due to reduced land cover. Accordingly, palynological proxies from S-Anatolia and S-Greece document a high human impact on the vegetation cover (Figure 5f-k). Thus, human activity became a dominant factor triggering geomorphodynamics. In addition to the general trend of climatic sensitivty to erosion, hydro-climatic signal from the Balkans and the area of modern Turkey differ. Whereas in the area of modern Turkey pronounced dry phases are recorded, an increase in wetness in the neighboring Balkans is recorded for the period between 3.0 ka BP and ca. 1.2 ka BP (Figure 5m-o).

Early Bronze Age to Iron Age
The reconstructed geomorphodynamics in the western lower Bakırçay plain decreased at the beginning of the Late Holocene. This is especially evident in the decline of the modeled sedimentation rate after 4.3 ka BP (Figure 5u). In addition, the cumulative probability of change-ages point to marked sediment facies change during that period (Figure 3e). The marked change in geomorphodynamics at the transition from the Middle to the Late Holocene is distinctly related to climatic variations: The so-called 4.2-ka-BP-drought-event caused a decline in settlement activities in the eastern Mediterranean (Figure 5a,e) and elsewhere around the Mediterranean Sea [117][118][119][120][121]); from Central Anatolia, it is well documented that many sites were abandoned [122]. The collapse of Mediterranean civilizations might also been fostered by the drought event [115,116,123].
In the western lower Bakırçay plain, an 800-year gap in occupation is recorded as from the end of the Early Bronze Age. During this period, the forest cover around Elaia on the coastal fringe of the western lower Bakırçay plain increased (Figure 5l).
The modeled sedimentation in the western lower Bakırçay plain remained on a relatively low level from the beginning of the Late Holocene (4.2 ka BP) until 3.3-3.2 ka BP. Although there is an increase in the cumulative probability of 14 C-ages directly after 4.2 ka BP, the observed CPF tendentiously remains on a low level until 2.7 ka BP (Figure 5r,s). After ca. 3 ka BP, geomorphodynamics in the western lower Bakırçay plain markedly increased, coinciding with a phase of low climatic erosion sensitivity and a general trend in increasing human activities that is also documented by an increase in palynological indicators for human activity ( Figure 5). The effect of settlement activities on geomorphodynamics is reflected by a relative increase in colluvial sediments after 3 ka BP (Figure 4). The period between 3.2 and 1.5 ka BP is known as the Beyşehir Occupation Phase in Anatolia. This phase was identified from a pollen diagram originating from southwestern Anatolia [107]. The meta-analysis of Woodbridge et al. [106] showed an increase in OJCV-pollen and grazing indicator species in Anatolia during the Beyşehir occupation phase. Pollen data indicate a contemporaneous decrease in forest cover around the Bay of Elaia (Figure 5l).
Not only in the western lower Bakırçay plain, but also in other regions in Anatolia, geomorphodynamics increased around 3 ka BP [8,114,124,125]. In the surroundings of the ancient city of Sagalassos in Southern Anatolia, geomorphodynamics increased simultaneously to a phase of increased deforestation and intensified cultivation and pastoralism [126,127].

Classic and Post-Classic Period(s)
Modeled sedimentation in the western lower Bakırçay plain peaks between 2.1 and 1.6 ka BP; the cumulative probability of 14 C-ages is the highest after 2.0 ka BP; change-ages point to facies change around the same time ( Figure 3). Thus, geomorphodynamics were high during the Classical period, reaching a maximum during Roman Imperial times ( Figure 6).
From an archaeological point of view, the transition from Hellenistic to Roman Imperial times is of special interest, as it is associated with a marked change in the settlement pattern in the western lower Bakırçay plain and a population increase in the major city of Pergamon. After the Eumenian expansion of Pergamon in the 2nd century BCE [128], the population density in the micro-region decreased, but the population in the cities of Pergamon and Elaia increased [128]. The city area of Pergamon quadrupled in the 2nd century BCE [129]. The area of the lower city further grew under the reign of Trajan (98)(99)(100)(101)(102)(103)(104)(105)(106)(107)(108)(109)(110)(111)(112)(113)(114)(115)(116)(117), including the location of residential quarters, major buildings, and bathes in the surroundings. In the 2nd century BCE, also the city of Elaia grew [128].
Although the sediment sequences were not obtained from the direct vicinity of the city of Pergamon, the Hellenistic-Roman transformation might have also affected remote areas in the micro-region. Historical reports (Galen,5,49) imply that around 180,000 people lived in the Pergamon Micro-Region in the 2nd century BCE. Population estimates based on historical architecture indicate a number of circa 34,000 people that lived in the city of Pergamon during that time [129]. Between Hellenistic and Roman Imperial times, the city area increased from circa 110 ha to circa 230 ha [129]. From approximations, it can be presumed that an important portion of the micro-region was under utilization to supply the urban and rural population of Pergamon and its micro-region in the 2nd century CE [14,130,131]. Although timber for the urban building program was felled. (Alternatively, large-scale imports would have been required.) Palynological data from Elaia indicate a relatively sparse forest cover during that period [41]. Although vegetation dynamics in the area of modern Turkey are well correlated with climate fluctuations [104], it is reasonable to assume that the high demand for food is a significant driver of forest cover reduction in the western lower Bakırçay plain. The dryness in the area of modern Turkey ( Figure 5)-and therefore a reduced vegetation cover-might have enhanced and accelerated the effect of human activities that lead to the increased geomorphodynamics since 2.0 ka BP.  Table A1); also the non-normalized data are displayed, as shown in Figure 3. Grey backgrounds indicate the variability within the archaeological period. The dashed lines indicated no difference between the observed data and the null model. SD = standard deviation, SR = sedimentation rate, N = Neolithic, CL = Chalcolithic, BA = Bronze Age, IA = Iron Age.
In accordance with our findings, data from sediments from the Bay of Elaia reported by Seeliger and co-authors indicate increased erosion during this time in the catchments around the city of Elaia (after 2.7-2.6 ka BP [41]; around 2.3-2.2 ka BP [37]). Schneider and co-authors did not, however, observe an increase in sediment dynamics during the Classical period analyzing single sediment sequences from the western lower Bakırçay plain [48,49]; they invoke terraces as a factor that minimized soil erosion processes during Roman Imperial times around the city of Atarneus [49]. We, however, observed a decline in geomorphodynamics after Roman Imperial times. Failures of terraces might thus not have caused increased erosion [132,133].
The increased geomorphodynamics during the Classic period appears more or less ubiquitously in the Aegean region [7,30,81,114,[134][135][136][137]. By contrast, the reconstructed sediment dynamics in the surroundings of Sagalassos [8] peaked during the Iron Age, before the major settlement phase in the Classic period and the related intensive land use [138]. Thus, the peak in sediment dynamics around Sagalassos appeared during the increase in human pressure on the landscape, but not during the heydays of the city when human pressure was presumably on its maximum. By contrast, sediment dynamics in the western lower Bakırçay plain, rather increased with increasing human activities. Several aspects may explain the differences in the development of the micro-regions of Sagalassos and Pergamon, including (i) the general variability in the record; (ii) the difference in location of the sediment sequences concerning sediment sinks or sources (sediments from valley fills vs. colluvial/fan sediments; vicinity to the coast in the western lower Bakırçay plain); (iii) soil conservation measures such as terracing [49]), or (iv) variability in landscape sensitivity to erosion, resilience, or trigger-reaction delays. Dusar et al. [8] analyze different colluvial and fluvial archives, whereas the sediments in our data set are mainly obtained from colluvial deposits or low order catchments. Archives from colluvial deposits or low order catchments reflect human disturbance more directly than alluvial deposits of larger drainage systems [8,10,139,140]. The signal from sediment archives obtained from the floodplain of higher order catchments-thus topographically lower compartments of the sediment cascade-are more likely affected by a time lack. The resilience time in e.g., colluvial archives is therefore much lower than in floodplain sediments.
Soil depletion [16] and a decreased erodibility due to higher stone coverage of soils after a period of erosion during the Iron Age around Sagalassos are a reason why geomorphodynamics might peak before the climax in human activities during the Classic period [12,141]. The depletion of soils on slopes that were sensitive to erosion might have caused a shift of agricultural activities and settlement locations to plain areas around Sagalassos that were less susceptible to erosion [12]. Preliminary data of Knitter and Ludwig [142] on the archaeological site distribution pattern in the Pergamon Micro-Region indicate a change of the preferred locations in the Classic period. While the preferred settlement location before Roman Imperial times were on summits, they moved to slopes during Roman Imperial times. Based on the available data, it can be inferred that Roman Imperial sites were located in places characterized by higher soil erosion sensitivity compared to older sites. This might also explain the increase in sediment dynamics recorded for the western lower Bakırçay plain from Hellenistic to Imperial times.
After Roman Imperial times, geomorphodynamics decreased in the western lower Bakırçay plain. Changes from the dominance of arable farming to animal husbandry around Elaia [41] might have been local. Rather, the decline in geomorphodynamics is attributed to an overall decline in human activities in the western lower Bakırçay plain.
The geomorphodynamics in the post-Classic periods are not further discussed in the paper, as archaeological research in the area focused on Prehistoric to Classic periods.

Conclusions
Data from sediment sequences obtained from the western lower Bakırçay plain were examined concerning geomorphodynamics mainly based on cumulative probability functions of 14 C-ages. We, therefore, modified current approaches to cope with certain biases in the data by comparing observed cumulative probabilities with simulated probabilities and by estimating sediment accumulation.
Both, cumulative probability functions and modeled sediment accumulation points to several phases of increased geomorphodynamics that were not identified in the analysis of single sequences. Inference statistics show that during these phases geomorphodynamics are statistically significantly different from randomly distributed data. Thus, although the number of samples in our data set (76 terrestrial ages) is limited compared to supra-regional studies, we obtained reasonable results that help to further understand human-environment interactions in the western lower Bakırçay plain.
Several aspects can be highlighted: • Phases of increased or reduced geomorphodynamics in the western lower Bakırçay plain follow the general trend of sediment dynamics in the eastern Mediterranean.
• The development of local geomorphodynamics is in good agreement with the changing hydro-climatic conditions as well as vegetation and population dynamics in Anatolia, Greece and the neighboring Balkan region. • From this coincidence, climatic triggers of geomorphodynamics appear to be (most) important in the Early Holocene, prior to Holocene aridization. However, the Early Holocene phase of increased geomorphodynamics coincided with the potential onset of early settlement activities in the western lower Bakırçay plain. Thus, relatively low human impact may have (even marginally) contributed to increased geomorphodynamics in the area-conceivably due to a general high sensitivity of the landscape to erosion. The decrease in geomorphodynamics in the western lower Bakırçay plain at the end of the middle Holocene is related to a climatic event, i.e., the 4.2-ka-BP-drought-event. This drought event caused widespread settlement abandonment and a decrease in population in Anatolia. • Geomorphodynamics in the western lower Bakırçay plain peaked during the Classic period. The transformation of the Pergamon Micro-Region from Hellenistic to Roman Imperial times included urbanization and demographic growth, which most likely were the most important triggers of increased geomorphodynamics. • Geomorphodynamics in the western lower Bakırçay plain and many other areas in Anatolia show a similar trend. Nonetheless, major differences in the Classical periods occurred. Whereas geomorphodynamics and settlement activities in the western lower Bakırçay plain are entangled, this is not the case in other major cities in Anatolia, such as Sagalassos. The reasons for the difference might either be related to differences between the sediment archives or diverging landscape sensitivities and land use dynamics.
To better understand the difference in geomorphodynamics and to overcome the simplicity of interpretations mainly based on synchronicity of geomorphodynamics and general trends in human activities and climate variability, further research on the entanglement of local geomorphodynamics and settlement dynamics is required. This especially includes the understanding of geomorphodynamics in the Early Holocene and studies on the cause-effect relationships between settlement activity and varying geomorphodynamics in the western lower Bakırçay plain [143]. This should include the modeling of past soil erosion to overcome the bias of location of the sediment sequences in the vicinity of settlements of a distinct period and the general problem of the spotty nature of sediment sequences. Evaluating the relationship between archaeological site patterns and landscape sensitivity to erosion, soil conservation measures, and studies that are not related to archaeological sites of a specific period might provide further insights in human-environment interactions in the western lower Bakırçay plain. In addition, palynological or phytolith studies might improve the understanding of agricultural use in the area.
Irrespective of the potential shortcomings of a meta-analysis, we conclude that our study reveals the general interrelationship of settlement history and geomorphodynamics in the western lower Bakırçay plain, especially the human-inducted climax of geomorphodynamcis in Roman Imperial times.

Acknowledgments:
Research was conducted within the project Die Transformation der Mikroregion Pergamon zwischen Hellenismus und römischer Kaiserzeit (The Transformation of the Micro-Region Pergamon between the Hellenistic Period and the Roman Imperial Age). Thanks to Xun Yang for commenting on a draft of the manuscript. We are grateful to three anonymous reviewers for their detailed comments that clearly helped to improve the article. We acknowledge support by the Open Access Publication Fund of the Freie Universität Berlin.

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.

Abbreviations
The following abbreviations are used in this manuscript: where SR i is the sedimentation rate related to a 14 C-age i, d i is the sampling depth of the age, and t i is the median date of the dating likelihood of the age. Taking the sedimentation rate, the maximum age of a sediment sequences is estimated (if only one date is available for a sequence; n i,j = 1):

BA
where t d max,j is the maximum age of the sediment sequence j and d max,j is the thickness of sediment sequence j. For sediment sequences with more than one 14 C-age available, we used a power function to estimate the maximum age: where e j is the power calculated by a nonlinear least squares method applied to all 14 C-ages available from a sediment sequence. The general form of the equation is: