- Email: [email protected]
The timing of the Black Sea flood event: Insights from modeling of glacial isostatic adjustment
Earth and Planetary Science Letters 452 (2016) 178–184
Contents lists available at ScienceDirect
Earth and Planetary Science Letters www.elsevier.com/locate/epsl
The timing of the Black Sea flood event: Insights from modeling of glacial isostatic adjustment Samuel L. Goldberg a,∗ , Harriet C.P. Lau a , Jerry X. Mitrovica a , Konstantin Latychev b a b
Department of Earth and Planetary Sciences, Harvard University, 20 Oxford Street, Cambridge, MA 02138, USA Department of Physics, University of Toronto, 60 St. George Street, Toronto, ON M5S 1A7, Canada
a r t i c l e
i n f o
Article history: Received 19 March 2016 Received in revised form 20 May 2016 Accepted 9 June 2016 Available online 16 August 2016 Editor: A. Yin Keywords: Black Sea sea level tectonics ice age Aegean Sea archaeology
a b s t r a c t We present a suite of gravitationally self-consistent predictions of sea-level change since Last Glacial Maximum (LGM) in the vicinity of the Bosphorus and Dardanelles straits that combine signals associated with glacial isostatic adjustment (GIA) and the flooding of the Black Sea. Our predictions are tuned to fit a relative sea level (RSL) record at the island of Samothrace in the north Aegean Sea and they include realistic 3-D variations in viscoelastic structure, including lateral variations in mantle viscosity and the elastic thickness of the lithosphere, as well as weak plate boundary zones. We demonstrate that 3-D Earth structure and the magnitude of the flood event (which depends on the pre-flood level of the lake) both have significant impact on the predicted RSL change at the location of the Bosphorus sill, and therefore on the inferred timing of the marine incursion. We summarize our results in a plot showing the predicted RSL change at the Bosphorus sill as a function of the timing of the flood event for different flood magnitudes up to 100 m. These results suggest, for example, that a flood event at 9 ka implies that the elevation of the sill was lowered through erosion by ∼14–21 m during, and after, the flood. In contrast, a flood event at 7 ka suggests erosion of ∼24–31 m at the sill since the flood. More generally, our results will be useful for future research aimed at constraining the details of this controversial, and widely debated geological event. © 2016 Elsevier B.V. All rights reserved.
1. Introduction The shallow Bosphorus and Dardanelles straits in northwest Turkey separate the Black Sea from the Mediterranean Sea and the global ocean system (Fig. 1). At the time of the Last Glacial Maximum (LGM), local sea level was low enough that the connection of the Black Sea to the ocean was broken (Clark et al., 2009). At this time, a freshwater or brackish lake, fed by riverine and glacial meltwater inputs, occupied the Black Sea basin (YankoHombach et al., 2014). At some time during post-LGM period, the rising waters of the Mediterranean breached the shallowest of the sills between the Mediterranean and Black Seas and the present connection between the two water bodies was established (Ryan et al., 2003). The nature of the reconnection remains uncertain. It may have been characterized by a gradual influx of seawater into the lake (Aksu et al., 2002; Hiscott and Aksu, 2002), by a catastrophic flood (Ryan et al., 1997, 2003; Ryan and Pitman, 1998; Major et al., 2002), or by an event of intermediate energy. In any case, the possibility that a cataclysmic natural event occurred has
*
Corresponding author. E-mail address: [email protected] (S.L. Goldberg).
http://dx.doi.org/10.1016/j.epsl.2016.06.016 0012-821X/© 2016 Elsevier B.V. All rights reserved.
potentially profound archaeological and historical implications for the Neolithic humans living in the area, and has been suggested to be the inspiration for the biblical flood legend (Ryan and Pitman, 1998). Several other aspects of the Black Sea flood event have also been sources of long-standing debate and controversy. The timing of the event, the location and height of the last sill to be breached, and the pre-flood level of the lake, are all uncertain. We review each of these issues, in turn. The shallowest of the presently existing sills between the Mediterranean and Black Seas lies at the southern entrance of the Bosphorus strait (see Fig. 1 and Fig. 2). The modern sill is 32–34 m below sea level, and consists of Quaternary sand overlying Paleozoic bedrock (Algan et al., 2001). In contrast, the sill depth at the Dardanelles strait is ∼80 m below sea level. The topography of the bedrock in the Bosphorus Strait features three sills 80–85 m below sea level. Sedimentation on these sills initiated prior to 10 ka and continued until 5.3 ka (Algan et al., 2001). For example, Holocene shells have been found at a depth of 70 m in sediment cores (Major et al., 2002). Thus, the formation of the modern sill must have continued through the period of post-glacial sea-level rise. Accordingly, there is no compelling evidence to sug-
S.L. Goldberg et al. / Earth and Planetary Science Letters 452 (2016) 178–184
179
Fig. 1. Topographic map of the eastern Mediterranean Sea and Black Sea region. Thin red lines indicate major plate boundary faults. The three blue dots indicate locations mentioned in the paper: (1) Bosphorus Strait, (2) Dardanelles Strait, and (3) Samothrace. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
Fig. 2. Schematic cross-section of the modern Black Sea sills, with respect to present-day sea level. The depths of the two sills are drawn to scale relative to each other; the depths of the basins in between are not. The Sea of Marmara has an average depth of ∼500 m; the Black Sea quickly drops to ∼2000 m depth beyond the Bosphorus.
gest that the elevations of either the modern sill or the bedrock sills is representative of the sill height at the time of the initial ocean incursion into the Black Sea. The timing of the flooding event is similarly poorly constrained. Soulet et al. (2011) analyzed microfossils from Black Sea sediments and estimated the initial appearance of marine species to be 9 ka. They adopted this value as an estimate of the timing of the sill breach. Yanko-Hombach et al. (2014) used several different paleosalinity proxies in sediment cores from the Black Sea, including macrofossils, microfossils, and pollen, to argue for a gradual marine transgression beginning no later than 8.9 ka. However, other studies have proposed a significantly earlier flooding event. Major et al. (2002) measured the isotopic compositions of mollusk shells, and on this basis suggested that the first marine influx took place at 12.8 ka. Aksu et al. (2002) examined the stratigraphy of the shelf sediments, using cores and seismic profiles, and concluded that the water level in the Black Sea Basin increased at 11–10.5 ka. Ryan et al. (1997) argued for a late, rapid flood at 7.1 ka based on dated shells from sediment cores. An additional uncertainty is the height of the water surface of the pre-flood lake in the Black Sea basin. Widely diverging estimates of this height, based on various geomorphological and sedimentological arguments, have appeared in the literature. For example, some have argued for a pre-flood surface ∼100 m below present sea level (Ryan, 2007; Lericolais et al., 2010; Nicholas et al., 2011), while others have concluded that the water surface was only 30–40 m below present sea level (Giosan et al., 2009; YankoHombach et al., 2014). It has also been suggested that the pre-
flood lake was not endorheic, but rather had sufficient inflow to crest the sill and flow into the lower Mediterranean Sea. This implies that the lake surface was at or just above the height of the sill (Algan et al., 2001) and that the reconnection during post-glacial sea level rise did not involve a rapid flood. Finally, there is also evidence that during periods of separation from the global ocean, the Black Sea lake alternated between a high mode with outflow into the Mediterranean during cold periods and an endorheic low mode with significant evaporation during warm periods (Ryan et al., 2003). Despite broad interest in the Black Sea flood hypothesis, only a single study by Lambeck et al. (2007) has presented geophysical predictions of ice age sea-level change (or glacial isostatic adjustment, GIA) within the region. (An additional study has considered the perturbation to Earth rotation driven by the flood event Spada et al., 1999.) The Lambeck et al. (2007) analysis computed sealevel changes across the last glacial cycle using an ice age model that was tuned to fit relative sea level (RSL) records at several tectonically stable Mediterranean sites, including: Carmel Coast, Israel; Villefranche, France; Versilia Plain, Italy; and Tiryns, Greece, in the Peloponnese. Their prediction of the sea-level history at the Bosphorus and Dardanelles sills allowed them to estimate the timing of the flood event for several different choices for the highest barrier at the time of the breach. For example, adopting the top of the present-day Bosphorus sill (currently at 32–34 m below sea level) as the highest barrier (i.e. assuming that the sill height was not modified by erosion or sedimentation during or after the flooding) they predicted that the connection of the Black Sea and the
180
S.L. Goldberg et al. / Earth and Planetary Science Letters 452 (2016) 178–184
Mediterranean took place between 9.5–10.3 ka. In contrast, using the top of the present-day depth of the Dardanelles sill (∼80 m) as the highest barrier yielded a predicted timing of ∼15 ka for the flood event. The sea-level predictions of Lambeck et al. (2007) have been adopted in subsequent studies to consider other flood scenarios. For example, Soulet et al. (2011) argued that the breach occurred at 9 ka at the Bosphorus sill, and using the Lambeck et al. (2007) results they estimated that the top of the sill height at this time was ∼10 m higher than the modern sill (assuming a present sill height of −35 m, and thus if this earlier sill height was preserved it would now be ∼25 m below sea level). The Lambeck et al. (2007) analysis represents an important initial study of ice age sea-level change in the region and the implications of this change for the timing of the connection between the Black Sea and the Mediterranean, but it had several shortcomings. First, the study did not consider the impact on the sea level predictions of the uncertainty in the pre-flood level of the Black Sea or in the timing of the flood. In particular, it assumed that water was continuously added to the Black Sea basin up to the level of the evolving local (computed) gravitational potential defining the sea-surface, effectively assuming that the Black Sea was always connected to the Mediterranean and the global ocean. This implies that the LGM level of the Black Sea was ∼130 m lower than present-day (or ∼100 m lower than the elevation of the present-day Bosphorus sill) and that at the time of the breach the sea surface on both sides of the Bosphorus sill was the same. Second, the GIA model adopted by Lambeck et al. (2007) assumed that Earth structure, including the thickness of the elastic lithosphere and mantle viscosity, varied with depth only. Since the region is a tectonically active plate boundary zone (Fig. 1), crust and mantle structure will have significant 3-D variability. Finally, the GIA model used in the analysis was tuned to fit Holocene sea level histories at sites that were relatively distant from the Bosphorus and Dardanelles straits. The closest such site (Tiryns, Greece) is located ∼430 km and ∼670 km from the Dardanelles and Bosphorus straits, respectively. In this study, we present a new sequence of predictions of GIA in the region based on the sea-level theory described in Kendall et al. (2004). Our analysis extends the analysis in Lambeck et al. (2007) in several important ways. First, we tune our ice age model to fit a new RSL history at a site in the northeast Aegean Sea (Vacchi et al., 2014). This site, the island of Samothrace, is 80 km from the Dardanelles strait and 310 km from the Bosphorus. Second, our GIA predictions incorporate three-dimensional crust and mantle structure using a finite-volume numerical approach described in Latychev et al. (2005). The geometry of this structure is constrained using a recently published, high-resolution seismic tomography model for the European upper mantle (Fichtner et al., 2013). Finally, our calculations will incorporate a range of pre-flood lake levels and sill depths. Our goal is to present revised predictions of the timing of the breach of the sill as a function of these two parameters. 2. Methods We use the gravitationally self-consistent sea-level theory of Kendall et al. (2004), which incorporates time-varying ocean geometry due to changes in the perimeter of grounded marine-based ice in polar regions and migration of shorelines due to sea-level changes. The theory requires, as inputs, models for both the global ice history across the last glacial cycle and the viscoelastic structure of the Earth’s mantle. In the case of the former, we adopt a form of the ICE-5G model of the last deglaciation phase (Peltier, 2004) that has been modified to fit a relative sea-level history near our region of interest, as explained below. We have not included
the feedback of Earth rotation into sea level, since this signal is small (∼1 m), and does not vary across the region of interest. Nearly all previous studies of GIA have adopted Maxwell viscoelastic Earth models with structure that varies with depth alone (e.g. Peltier, 2004; Lambeck et al., 1998; Mitrovica and Forte, 2004). The Bosphorus and Dardanelles straits are located close to major plate boundaries (Fig. 1) and this suggests that both the crustlithospheric thickness and mantle viscosity beneath the region will exhibit significant lateral variations. Accordingly, our GIA calculations are based on a finite-volume numerical methodology that incorporates large amplitude, three-dimensional lateral variability in the Earth’s rheological parameters (Latychev et al., 2005). The numerical scheme discretizes the Earth into a global tetrahedral grid, and recent improvements in the software allow for regional grid refinement. We have taken advantage of the latter, and within the Mediterranean and Black Seas the grid has a surface spatial resolution of 6 km. Outside this region the resolution drops to 12 km. Lateral variations in elastic lithospheric thickness of the Earth model are prescribed from the model of Conrad and Lithgow-Bertelloni (2006), and within this variability plate boundaries are incorporated as a 100 km zone of low viscosity to allow decoupling of adjacent plates. The mantle viscosity field is constructed by mapping tomographic models of seismic shear wave velocity heterogeneity into lateral viscosity variations using the methodology described in Latychev et al. (2005). For this purpose, we adopt the seismic model S40RTS (Ritsema et al., 2011) to first build a global mantle viscosity field, and then replace the upper mantle viscosity below Europe with structure mapped from a highresolution tomographic model of seismic structure in the region (Fichtner et al., 2013). The result is an Earth model with a ∼3 orders of magnitude change in viscosity beneath the region encompassing the western Black Sea and eastern Aegean Sea (Fig. 3). The spherically averaged value of the 3-D Earth model is characterized by a lithospheric thickness of ∼100 km, and upper and lower mantle viscosities of 5 × 1020 Pa s and 5 × 1021 Pa s, respectively. 3. Results GIA will significantly affect the post-glacial sea-level signal in the region covering the Bosphorus and Dardanelles straits and our prediction of this impact should be consistent with geological records of relative sea-level (RSL) change in the vicinity of this area. As we have discussed, Lambeck et al. (2007) tuned their GIA simulation to be consistent with RSL records at four sites, the closest of which was 430–670 km from the two sills. Subsequent to their work, Vacchi et al. (2014) published RSL data from the northeast Aegean Sea and we have revised the ice history in the GIA simulation in order to fit the record from the site closest to the Bosphorus and Dardanelles straits, the island of Samothrace (Fig. 1). The sea level data were obtained from a site on the western tip of the island (Vacchi et al., 2014). Specifically, the ICE-5G ice history was modified by shallowing the eustatic curve by ∼2 m at 8 ka, and tapering this perturbation to 0 m by 4 ka, in order to provide a better fit to the local data set. This revised ice history, in combination with the 3-D Earth model described above, yields a good fit to the sea-level record at Samothrace (Fig. 4, red lines). As a point of comparison, we also show in Fig. 4 an RSL prediction at Samothrace based on the same ice history but using an Earth model that varies with depth alone (blue line). The model is identical to the spherically averaged rheology of the 3-D Earth model. The prediction of this model exhibits a much shallower sea-level trend since 10 ka, and significantly misfits the geological record of sea-level change. We performed a sequence of simulations based on a large suite of 1-D Earth models in which we varied the lithospheric thickness, upper mantle viscosity and lower mantle viscosity and found that none of the predictions using 1-D
S.L. Goldberg et al. / Earth and Planetary Science Letters 452 (2016) 178–184
181
Fig. 3. Lateral variations in mantle viscosity at a depth of 200 km below the western Black Sea and Aegean Sea incorporated into the GIA predictions described in the text. The symbol ν3D denotes the viscosity adopted in the calculation, and ν1D is the spherically averaged value of viscosity at this depth (5 × 1020 Pa s). The inset displays the variation in lithospheric thickness across the same region, including the location of plate boundaries (blue lines), incorporated into the lithospheric model. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
Fig. 4. Predictions of relative sea level at Samothrace (Fig. 1) for different GIA simulations. Red line: RSL prediction based on the ice history and 3-D Earth model discussed in the text, but no flooding of the Black Sea. Dashed, dotted red lines: Same as solid red line, except that flooding of the Black Sea that increased its bathymetry by an average of 50 m and 100 m, respectively, is included in the calculation. Blue line: Same as solid red line, except that lateral variations in Earth structure are ignored. The observations all involve the marine bivalve Cerastoderma glaucum and are taken from Vacchi et al. (2014). Two older observations reported in that study are not included. These are based on marine organics, represent terrestrial limits only and have potentially large age errors (A. Rovere, personal communication, 2016). The RMS misfit between the 3-D calculations and the observations is a factor of ∼3.2 smaller than the misfit associated with the 1-D calculation. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
Earth models provided a reasonable fit to the observations (see supplement, Fig. S1), suggesting that incorporating lateral variations in Earth structure is necessary to fit the observed data. To explore this issue further, Fig. 5 shows the difference in the RSL predictions at 10 ka over southern Europe and north Africa between the 3-D and 1-D Earth models. These differences are significant, with peak amplitudes of ∼10 m. The 3-D prediction is clearly influenced by large-scale variability in mantle viscosity and
lithospheric thickness; however, the impact of plate boundary effects is also noticeable in localized structure that extends across the region. For example, the southern boundary of the Eurasian plate is evident in a localized perturbation in sea level that extends west across Greece, and along the northern boundaries of the Aegean and Anatolian Plates, passing in the vicinity of the Dardanelles and Bosphorus straits. Since the boundary is modeled as a low viscosity zone within the lithosphere, the 3-D GIA calculation permits the crust on either side of the boundary to move relatively independently in response to the water loading that extends across the Mediterranean and Aegean. This permits greater subsidence of the Aegean crust (and thus sea level rise) in response to post-LGM meltwater loading, relative to the continuous lithosphere of the 1-D earth model prediction, which explains the results in Fig. 4. Two other sea level records exist in the Aegean Sea, further from the sills, on the islands of Limnos and Iskele (Vacchi et al., 2014). The fit of the 3-D model to these curves is less satisfactory than at Samothrace (see supplement, Figure S2). However, the incorporation of 3-D structure progressively improves the fit relative to the 1-D model as one considers sites closer to the sill. Given the short-wavelength variability evident in Fig. 5, we have chosen to tune the model to fit the site closest to the Black Sea sills. The solid red and blue lines in Fig. 4 and Fig. S2 were computed with the Black Sea masked so that no meltwater was allowed to enter it during the deglaciation phase. This raises the issue as to whether the flooding of the Black Sea would impact the sealevel trend predicted at Samothrace, and whether our tuning of the regional GIA model should account for this. To consider this, we ran simulations in which we flooded the Black Sea with varying volumes of water and tracked the resulting regional sea-level change over the next 10 kyr. As an example, Fig. 6 shows the total sea-level change 10 kyr after a flooding event that increased the bathymetry of Black Sea level by an average of 100 m. While the predicted impact of the flood on sea level is significant at the Bosphorus strait and sill, it is relatively minor at the Dardanelles sill (∼1 m) and at Samothrace (∼0.5 m). This is verified in Fig. 4, where the dashed and dotted red lines represent the predicted RSL curves at Samothrace for scenarios in which the original GIA calculation (solid red line) was augmented to include a Black Sea flood
182
S.L. Goldberg et al. / Earth and Planetary Science Letters 452 (2016) 178–184
Fig. 5. Difference in the predicted RSL at 10 ka (m) due to GIA based on calculations that adopt the 3-D and 1-D viscoelastic Earth models discussed in the text (i.e., the 3-D minus 1-D RSL predictions).
Fig. 6. Predicted total change in sea level (m) 10 kyr after a flooding event that raised the average bathymetry of the Black Sea by 100 m. The color on the main figure saturates at 12 m; the maximum value is ∼16 m in the center of the Black Sea. The prediction adopts the 3-D viscoelastic Earth model described in the text.
at 10 ka equivalent to an average bathymetry change of 50 m and 100 m, respectively. We conclude that our tuning of the GIA model to fit the RSL record at Samothrace is not sensitive to the timing or magnitude of the flooding event. However, in contrast to the predictions at Samothrace, the prediction in Fig. 6 indicates that the flooding event had a significant impact on the sea level change at the Bosphorus sill, and must therefore be taken into account when estimating the timing or magnitude of the flood event. To quantify this sensitivity, Fig. 7 shows predictions of GIA-induced RSL at Bosphorus based on the ice history discussed above and including various scenarios for the magnitude (Fig. 7A) and timing (Fig. 7B) of the Black Sea flood
event. The scenarios in Fig. 7A span the range of flood magnitudes proposed in the literature, and it is clear that these predictions are sensitive to this uncertainty. Indeed, relative to the no flood case, the 100 m event considered in Fig. 7A leads to a ∼10 m greater rise in the sea-level change predicted at the Bosphorus sill from 10 ka to present. This increased sea-level rise is largely due to crustal subsidence at the sill driven by the incursion into the Black Sea and associated water loading. Our complete set of calculations is used to construct Fig. 8, which plots RSL change at the location of the sill since the flood event versus the time of the event for several different flood magnitudes (in units of mean change in Black Sea depth). This figure
S.L. Goldberg et al. / Earth and Planetary Science Letters 452 (2016) 178–184
183
Fig. 7. Predictions of RSL change at the Bosphorus sill due to GIA and the Black Sea flood event. (A) Computed RSL change in the case of a flood event at 10 ka of varying magnitude (as labeled). (B) Computed RSL for a flood event of variable timing (as labeled) that raised the average depth of the Black Sea by 100 m. All calculations adopt the 3-D viscoelastic Earth model discussed in the text (see Fig. 3).
The results in Fig. 8 allow us to reassess previously published scenarios for the Black Sea flood event. For example, based on a synthesis of micropaleontological reconstructions, Soulet et al. (2011) estimated that the flood event occurred at 9 ka and, using GIA results in Lambeck et al. (2007), concluded that the sill has experienced 10 m of erosion during, and subsequent to, the flooding. Our results suggest that for a flood at 9 ka, erosion of 14–21 m occurred to bring the sill height to its present level of −33 m. If instead we adopt the 7.1 ka timing of Ryan et al. (1997), we estimate 24–30 m of erosion of the sill (i.e. if the sill had not eroded, it would be close to present sea level). In both cases, the range in the estimate reflects the uncertainty in the magnitude of the flood. 4. Discussion and conclusions
Fig. 8. Predicted RSL at the location of the Bosphorus sill at the time of imposed flooding for three different cases of flood magnitude (as labeled). Predictions are based on the ice history and 3-D Earth model described in the text.
allows for the estimation of one of the flood parameters (timing, magnitude, or sill depth) as a function of the other two. Consider, as an example, a flood event at 10 ka of magnitude 100 m. In this case, the RSL change since the flood is predicted to be −25 m. This prediction changes to −20.5 m and −17 m for flood events with the same timing and of magnitude 60 m and 30 m, respectively. Since the present day depth of the sill is ∼33 m, these scenarios imply that erosion of the sill over the last 10 kyr has lowered its elevation by 8 m, 12.5 m, and 16 m, for flood magnitudes of 100 m, 60 m, and 30 m, respectively. More generally, any scenario in Fig. 8 that yields a RSL above −33 m implies that the sill has eroded down to the observed level of −33 m since the time of flooding, while any scenario with a predicted RSL below −33 m implies that sedimentation has raised the elevation of the sill over the same time period. Thus, the results in Fig. 8 indicate that the sill must have experienced erosion if the flood event occurred at, or after, ∼11 ka. If instead the sill depth has not changed since the flood event, the event occurred between 11 and 11.6 ka.
We have predicted relative sea-level change at the location of the Bosphorus sill and surrounding regions since LGM using a gravitationally self-consistent ice age sea-level theory and a GIA model tuned to fit a RSL record in the northern Aegean. Our results indicate that the numerical sea-level predictions are sensitive to the incorporation of 3-D viscoelastic Earth structure, including lateral variations in mantle viscosity, lithospheric thickness and plate boundaries, as well as to the magnitude and timing of the flooding of the Black Sea. Neither of these sensitivities, which can perturb the sea-level predictions by as much as ∼10 m, was considered in the only previous GIA modeling of sea-level change in the region. We note that our interpretation of the RSL predictions has assumed tectonic stability over the time interval of study. While the area surrounding the Bosphorus is tectonically active, the magnitude of vertical tectonic motion that may have occurred over the past 10 ka is thought to be small. In the western Marmara Sea, which is a more tectonically active area, late Quaternary uplift rates have been calculated at ∼0.3 mm/yr, generally decreasing eastwards towards the Bosphorus (Yaltirak et al., 2002). Close to the Bosphorus itself, late Miocene limestones have been uplifted to 200 m, which indicates an average uplift rate of ∼0.02 mm/yr (Gökasan et al., 1997). This rate suggests an uplift of ∼0.2 m since 10 ka assuming no recent acceleration of the uplift, which is well within the uncertainties of our predictions. We have summarized our predictions in a plot showing RSL change at the location of the Bosphorus sill as a function of the
184
S.L. Goldberg et al. / Earth and Planetary Science Letters 452 (2016) 178–184
timing of the Black Sea flood for a suite of different flood magnitudes ranging from 30–100 m. As we noted in the Introduction, a wide range of estimates exist for the timing of the Black Sea flood and the issue has been the subject of significant, ongoing controversy – a debate that no doubt reflects the possible connection between the flood event and the Biblical flood legend. Our goal in the present study has been to provide a state-of-the-art geophysical analysis of this event that explicitly accounts for existing uncertainties in the nature of the sill and lake elevation at the time of the possibly cataclysmic event. Acknowledgements We thank two anonymous reviewers for constructive comments on a previous version of the manuscript. This work was supported by a HUCE summer research grant (to SLG), Harvard University, and the Canadian Institute for Advanced Research. Appendix A. Supplementary material Supplementary material related to this article can be found online at http://dx.doi.org/10.1016/j.epsl.2016.06.016. References Aksu, A.E., Hiscott, R.N., Yasar, D., Isler, F.I., Marsh, S., 2002. Seismic stratigraphy of Late Quaternary deposits from the southwestern Black Sea shelf: evidence for non-catastrophic variations in sea-level during the last ∼10000 years. Mar. Geol. 190, 61–94. Algan, O., Cagatay, N., Tchepalyga, A., Ongan, D., Eastoe, C., Gokasan, E., 2001. Stratigraphy of the sediment infill in Bosphorus Strait: water exchange between the Black and Mediterranean Seas during the last glacial Holocene. Geo Mar. Lett. 20, 209–218. Clark, P., Dyke, A., Shakun, J., Carlson, A., Clark, J., Wohlfarth, B., Mitrovica, J., Hostetler, S., McCabe, A.M., 2009. The Last Glacial Maximum. Science 325, 710–714. Conrad, C.P., Lithgow-Bertelloni, C., 2006. Influence of continental roots and asthenosphere on plate–mantle coupling. Geophys. Res. Lett. 33, L05312. Fichtner, A., Saygin, E., Tayaz, T., Cupilard, P., Capdeville, Y., Trampert, J., 2013. The deep structure of the North Anatolian Fault Zone. Earth Planet. Sci. Lett. 373, 109–117. Giosan, L., Filip, F., Constatinescu, S., 2009. Was the Black Sea catastrophically flooded in the early Holocene? Quat. Sci. Rev. 28, 1–6. Gökasan, E., Demirbag, E., Oktay, F., Ecevitoglu, B., Simsek, M., Yuce, H., 1997. On the origin of the Bosphorus. Mar. Geol. 140, 183–199. Hiscott, R.N., Aksu, A.E., 2002. Late Quaternary history of the Marmara Sea and Black Sea from high-resolution seismic and gravity-core studies. Mar. Geol. 190, 261–282. Kendall, R.A., Mitrovica, J.X., Milne, G.A., 2004. On post-glacial sea level – II. Numerical formulation and comparative results on spherically symmetric models. Geophys. J. Int. 161, 679–706.
Lambeck, K., Sivan, D., Purcell, A., 2007. Timing of the last Mediterranean Sea–Black Sea connection from isostatic models and regional sea-level data. In: YankoHombach, V.V., Gilbert, A.S., Panin, N., Dolukhanov, P.M. (Eds.), The Black Sea Flood Question: Change in Coastline, Climate and Human Settlement. Springer, Heidelberg, pp. 797–808. Lambeck, K., Smither, C., Johnston, P., 1998. Sea-level change, glacial rebound and mantle viscosity for northern Europe. Geophys. J. Int. 134, 102–144. Latychev, K., Mitrovica, J., Tromp, J., Tamisiea, M., Komatitsch, D., Christara, C., 2005. Glacial isostatic adjustment on 3-D Earth models: a finite-volume formulation. Geophys. J. Int. 161, 421–444. Lericolais, G., Guichard, F., Morigi, C., Minereau, A., Popescu, I., Radan, S., 2010. A post Younger Dryas Black Sea regression identified from sequence stratigraphy correlated to core analysis and dating. Quat. Int. 225 (2), 199–209. Major, C., Ryan, W., Lericolais, G., Hajdas, I., 2002. Constraints on Black Sea outflow to the Sea of Marmara during the last glacial–interglacial transition. Mar. Geol. 190, 19–34. Mitrovica, J., Forte, A., 2004. A new inference of mantle viscosity base upon joint inversion of convection and glacial isostatic adjustment data. Earth Planet. Sci. Lett. 225, 177–189. Nicholas, W.A., Chivas, A.R., Murray-Wallace, C.V., Fink, D., 2011. Prompt transgression and gradual salinisation of the Black Sea during the early Holocene constrained by amino acid racemization and radiocarbon dating. Quat. Sci. Rev. 30, 3769–3790. Peltier, W., 2004. Global glacial isostasy and the surface of the ice-age Earth: the ICE-5G (VM2) model and GRACE. Annu. Rev. Earth Planet. Sci. 32, 111–149. Ritsema, J., Deuss, A., van Heijst, H., Woodhouse, J., 2011. S40RTS: a degree-40 shear-velocity model for the mantle from new Rayleigh wave dispersion, teleseismic traveltime and normal-mode splitting function measurements. Geophys. J. Int. 184 (3), 1223–1236. Ryan, W.B.F., 2007. Status of the Black Sea flood hypothesis. In: Yanko-Hombach, V., Gilbert, A.S., Panin, N., Dolukhanov, P.M. (Eds.), The Black Sea Flood Question: Changes in Coastline, Climate, and Human Settlement. Springer, Dordrecht, pp. 3–88. Ryan, W.B.F., Pitman III, W.C., Major, C.O., Shimkus, K., Moskalenko, V., Jones, G.A., Dimitrov, P., Görür, N., Sakinç, M., Yüce, H., 1997. An abrupt drowning of the Black Sea shelf. Mar. Geol. 138, 119–126. Ryan, W., Pitman, W., 1998. Noah’s Flood: The New Scientific Discoveries About Events That Changed History. Simon and Schuster, New York. Ryan, W., Major, C., Lericolais, Gl., Goldstein, S., 2003. Catastrophic flooding of the Black Sea. Annu. Rev. Earth Planet. Sci. 31, 525–554. Soulet, G., Ménot, G., Lericolais, G., Bard, E., 2011. A revised calendar age for the last reconnection of the Black Sea to the global ocean. Quat. Sci. Rev. 30, 1019–1026. Spada, G., Alfonsi, L., Boschi, E., 1999. Chandler wobble excitation by catastrophic flooding of the Black Sea. Ann. Geofis. 42, 749–754. Vacchi, M., Rovere, A., Chatzipetros, A., Zouros, N., Firpo, M., 2014. An updated database of Holocene relative sea level changes in NE Aegean Sea. Quat. Int. 328–329, 301–310. Yaltirak, C., Sakinc, M., Aksu, A., Hiscott, R., Galleb, B., Ulgen, U., 2002. Late Pleistocene uplift history along the southwestern Marmara Sea determined from raised coastal deposits and global sea-level variations. Mar. Geol. 190, 283–305. Yanko-Hombach, V., Mudie, P., Kadurin, S., Larchenkov, E., 2014. Holocene marine transgression in the Black Sea: new evidence from the northwestern Black Sea shelf. Quat. Int. 345, 100–118.
Contents lists available at ScienceDirect
Earth and Planetary Science Letters www.elsevier.com/locate/epsl
The timing of the Black Sea flood event: Insights from modeling of glacial isostatic adjustment Samuel L. Goldberg a,∗ , Harriet C.P. Lau a , Jerry X. Mitrovica a , Konstantin Latychev b a b
Department of Earth and Planetary Sciences, Harvard University, 20 Oxford Street, Cambridge, MA 02138, USA Department of Physics, University of Toronto, 60 St. George Street, Toronto, ON M5S 1A7, Canada
a r t i c l e
i n f o
Article history: Received 19 March 2016 Received in revised form 20 May 2016 Accepted 9 June 2016 Available online 16 August 2016 Editor: A. Yin Keywords: Black Sea sea level tectonics ice age Aegean Sea archaeology
a b s t r a c t We present a suite of gravitationally self-consistent predictions of sea-level change since Last Glacial Maximum (LGM) in the vicinity of the Bosphorus and Dardanelles straits that combine signals associated with glacial isostatic adjustment (GIA) and the flooding of the Black Sea. Our predictions are tuned to fit a relative sea level (RSL) record at the island of Samothrace in the north Aegean Sea and they include realistic 3-D variations in viscoelastic structure, including lateral variations in mantle viscosity and the elastic thickness of the lithosphere, as well as weak plate boundary zones. We demonstrate that 3-D Earth structure and the magnitude of the flood event (which depends on the pre-flood level of the lake) both have significant impact on the predicted RSL change at the location of the Bosphorus sill, and therefore on the inferred timing of the marine incursion. We summarize our results in a plot showing the predicted RSL change at the Bosphorus sill as a function of the timing of the flood event for different flood magnitudes up to 100 m. These results suggest, for example, that a flood event at 9 ka implies that the elevation of the sill was lowered through erosion by ∼14–21 m during, and after, the flood. In contrast, a flood event at 7 ka suggests erosion of ∼24–31 m at the sill since the flood. More generally, our results will be useful for future research aimed at constraining the details of this controversial, and widely debated geological event. © 2016 Elsevier B.V. All rights reserved.
1. Introduction The shallow Bosphorus and Dardanelles straits in northwest Turkey separate the Black Sea from the Mediterranean Sea and the global ocean system (Fig. 1). At the time of the Last Glacial Maximum (LGM), local sea level was low enough that the connection of the Black Sea to the ocean was broken (Clark et al., 2009). At this time, a freshwater or brackish lake, fed by riverine and glacial meltwater inputs, occupied the Black Sea basin (YankoHombach et al., 2014). At some time during post-LGM period, the rising waters of the Mediterranean breached the shallowest of the sills between the Mediterranean and Black Seas and the present connection between the two water bodies was established (Ryan et al., 2003). The nature of the reconnection remains uncertain. It may have been characterized by a gradual influx of seawater into the lake (Aksu et al., 2002; Hiscott and Aksu, 2002), by a catastrophic flood (Ryan et al., 1997, 2003; Ryan and Pitman, 1998; Major et al., 2002), or by an event of intermediate energy. In any case, the possibility that a cataclysmic natural event occurred has
*
Corresponding author. E-mail address: [email protected] (S.L. Goldberg).
http://dx.doi.org/10.1016/j.epsl.2016.06.016 0012-821X/© 2016 Elsevier B.V. All rights reserved.
potentially profound archaeological and historical implications for the Neolithic humans living in the area, and has been suggested to be the inspiration for the biblical flood legend (Ryan and Pitman, 1998). Several other aspects of the Black Sea flood event have also been sources of long-standing debate and controversy. The timing of the event, the location and height of the last sill to be breached, and the pre-flood level of the lake, are all uncertain. We review each of these issues, in turn. The shallowest of the presently existing sills between the Mediterranean and Black Seas lies at the southern entrance of the Bosphorus strait (see Fig. 1 and Fig. 2). The modern sill is 32–34 m below sea level, and consists of Quaternary sand overlying Paleozoic bedrock (Algan et al., 2001). In contrast, the sill depth at the Dardanelles strait is ∼80 m below sea level. The topography of the bedrock in the Bosphorus Strait features three sills 80–85 m below sea level. Sedimentation on these sills initiated prior to 10 ka and continued until 5.3 ka (Algan et al., 2001). For example, Holocene shells have been found at a depth of 70 m in sediment cores (Major et al., 2002). Thus, the formation of the modern sill must have continued through the period of post-glacial sea-level rise. Accordingly, there is no compelling evidence to sug-
S.L. Goldberg et al. / Earth and Planetary Science Letters 452 (2016) 178–184
179
Fig. 1. Topographic map of the eastern Mediterranean Sea and Black Sea region. Thin red lines indicate major plate boundary faults. The three blue dots indicate locations mentioned in the paper: (1) Bosphorus Strait, (2) Dardanelles Strait, and (3) Samothrace. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
Fig. 2. Schematic cross-section of the modern Black Sea sills, with respect to present-day sea level. The depths of the two sills are drawn to scale relative to each other; the depths of the basins in between are not. The Sea of Marmara has an average depth of ∼500 m; the Black Sea quickly drops to ∼2000 m depth beyond the Bosphorus.
gest that the elevations of either the modern sill or the bedrock sills is representative of the sill height at the time of the initial ocean incursion into the Black Sea. The timing of the flooding event is similarly poorly constrained. Soulet et al. (2011) analyzed microfossils from Black Sea sediments and estimated the initial appearance of marine species to be 9 ka. They adopted this value as an estimate of the timing of the sill breach. Yanko-Hombach et al. (2014) used several different paleosalinity proxies in sediment cores from the Black Sea, including macrofossils, microfossils, and pollen, to argue for a gradual marine transgression beginning no later than 8.9 ka. However, other studies have proposed a significantly earlier flooding event. Major et al. (2002) measured the isotopic compositions of mollusk shells, and on this basis suggested that the first marine influx took place at 12.8 ka. Aksu et al. (2002) examined the stratigraphy of the shelf sediments, using cores and seismic profiles, and concluded that the water level in the Black Sea Basin increased at 11–10.5 ka. Ryan et al. (1997) argued for a late, rapid flood at 7.1 ka based on dated shells from sediment cores. An additional uncertainty is the height of the water surface of the pre-flood lake in the Black Sea basin. Widely diverging estimates of this height, based on various geomorphological and sedimentological arguments, have appeared in the literature. For example, some have argued for a pre-flood surface ∼100 m below present sea level (Ryan, 2007; Lericolais et al., 2010; Nicholas et al., 2011), while others have concluded that the water surface was only 30–40 m below present sea level (Giosan et al., 2009; YankoHombach et al., 2014). It has also been suggested that the pre-
flood lake was not endorheic, but rather had sufficient inflow to crest the sill and flow into the lower Mediterranean Sea. This implies that the lake surface was at or just above the height of the sill (Algan et al., 2001) and that the reconnection during post-glacial sea level rise did not involve a rapid flood. Finally, there is also evidence that during periods of separation from the global ocean, the Black Sea lake alternated between a high mode with outflow into the Mediterranean during cold periods and an endorheic low mode with significant evaporation during warm periods (Ryan et al., 2003). Despite broad interest in the Black Sea flood hypothesis, only a single study by Lambeck et al. (2007) has presented geophysical predictions of ice age sea-level change (or glacial isostatic adjustment, GIA) within the region. (An additional study has considered the perturbation to Earth rotation driven by the flood event Spada et al., 1999.) The Lambeck et al. (2007) analysis computed sealevel changes across the last glacial cycle using an ice age model that was tuned to fit relative sea level (RSL) records at several tectonically stable Mediterranean sites, including: Carmel Coast, Israel; Villefranche, France; Versilia Plain, Italy; and Tiryns, Greece, in the Peloponnese. Their prediction of the sea-level history at the Bosphorus and Dardanelles sills allowed them to estimate the timing of the flood event for several different choices for the highest barrier at the time of the breach. For example, adopting the top of the present-day Bosphorus sill (currently at 32–34 m below sea level) as the highest barrier (i.e. assuming that the sill height was not modified by erosion or sedimentation during or after the flooding) they predicted that the connection of the Black Sea and the
180
S.L. Goldberg et al. / Earth and Planetary Science Letters 452 (2016) 178–184
Mediterranean took place between 9.5–10.3 ka. In contrast, using the top of the present-day depth of the Dardanelles sill (∼80 m) as the highest barrier yielded a predicted timing of ∼15 ka for the flood event. The sea-level predictions of Lambeck et al. (2007) have been adopted in subsequent studies to consider other flood scenarios. For example, Soulet et al. (2011) argued that the breach occurred at 9 ka at the Bosphorus sill, and using the Lambeck et al. (2007) results they estimated that the top of the sill height at this time was ∼10 m higher than the modern sill (assuming a present sill height of −35 m, and thus if this earlier sill height was preserved it would now be ∼25 m below sea level). The Lambeck et al. (2007) analysis represents an important initial study of ice age sea-level change in the region and the implications of this change for the timing of the connection between the Black Sea and the Mediterranean, but it had several shortcomings. First, the study did not consider the impact on the sea level predictions of the uncertainty in the pre-flood level of the Black Sea or in the timing of the flood. In particular, it assumed that water was continuously added to the Black Sea basin up to the level of the evolving local (computed) gravitational potential defining the sea-surface, effectively assuming that the Black Sea was always connected to the Mediterranean and the global ocean. This implies that the LGM level of the Black Sea was ∼130 m lower than present-day (or ∼100 m lower than the elevation of the present-day Bosphorus sill) and that at the time of the breach the sea surface on both sides of the Bosphorus sill was the same. Second, the GIA model adopted by Lambeck et al. (2007) assumed that Earth structure, including the thickness of the elastic lithosphere and mantle viscosity, varied with depth only. Since the region is a tectonically active plate boundary zone (Fig. 1), crust and mantle structure will have significant 3-D variability. Finally, the GIA model used in the analysis was tuned to fit Holocene sea level histories at sites that were relatively distant from the Bosphorus and Dardanelles straits. The closest such site (Tiryns, Greece) is located ∼430 km and ∼670 km from the Dardanelles and Bosphorus straits, respectively. In this study, we present a new sequence of predictions of GIA in the region based on the sea-level theory described in Kendall et al. (2004). Our analysis extends the analysis in Lambeck et al. (2007) in several important ways. First, we tune our ice age model to fit a new RSL history at a site in the northeast Aegean Sea (Vacchi et al., 2014). This site, the island of Samothrace, is 80 km from the Dardanelles strait and 310 km from the Bosphorus. Second, our GIA predictions incorporate three-dimensional crust and mantle structure using a finite-volume numerical approach described in Latychev et al. (2005). The geometry of this structure is constrained using a recently published, high-resolution seismic tomography model for the European upper mantle (Fichtner et al., 2013). Finally, our calculations will incorporate a range of pre-flood lake levels and sill depths. Our goal is to present revised predictions of the timing of the breach of the sill as a function of these two parameters. 2. Methods We use the gravitationally self-consistent sea-level theory of Kendall et al. (2004), which incorporates time-varying ocean geometry due to changes in the perimeter of grounded marine-based ice in polar regions and migration of shorelines due to sea-level changes. The theory requires, as inputs, models for both the global ice history across the last glacial cycle and the viscoelastic structure of the Earth’s mantle. In the case of the former, we adopt a form of the ICE-5G model of the last deglaciation phase (Peltier, 2004) that has been modified to fit a relative sea-level history near our region of interest, as explained below. We have not included
the feedback of Earth rotation into sea level, since this signal is small (∼1 m), and does not vary across the region of interest. Nearly all previous studies of GIA have adopted Maxwell viscoelastic Earth models with structure that varies with depth alone (e.g. Peltier, 2004; Lambeck et al., 1998; Mitrovica and Forte, 2004). The Bosphorus and Dardanelles straits are located close to major plate boundaries (Fig. 1) and this suggests that both the crustlithospheric thickness and mantle viscosity beneath the region will exhibit significant lateral variations. Accordingly, our GIA calculations are based on a finite-volume numerical methodology that incorporates large amplitude, three-dimensional lateral variability in the Earth’s rheological parameters (Latychev et al., 2005). The numerical scheme discretizes the Earth into a global tetrahedral grid, and recent improvements in the software allow for regional grid refinement. We have taken advantage of the latter, and within the Mediterranean and Black Seas the grid has a surface spatial resolution of 6 km. Outside this region the resolution drops to 12 km. Lateral variations in elastic lithospheric thickness of the Earth model are prescribed from the model of Conrad and Lithgow-Bertelloni (2006), and within this variability plate boundaries are incorporated as a 100 km zone of low viscosity to allow decoupling of adjacent plates. The mantle viscosity field is constructed by mapping tomographic models of seismic shear wave velocity heterogeneity into lateral viscosity variations using the methodology described in Latychev et al. (2005). For this purpose, we adopt the seismic model S40RTS (Ritsema et al., 2011) to first build a global mantle viscosity field, and then replace the upper mantle viscosity below Europe with structure mapped from a highresolution tomographic model of seismic structure in the region (Fichtner et al., 2013). The result is an Earth model with a ∼3 orders of magnitude change in viscosity beneath the region encompassing the western Black Sea and eastern Aegean Sea (Fig. 3). The spherically averaged value of the 3-D Earth model is characterized by a lithospheric thickness of ∼100 km, and upper and lower mantle viscosities of 5 × 1020 Pa s and 5 × 1021 Pa s, respectively. 3. Results GIA will significantly affect the post-glacial sea-level signal in the region covering the Bosphorus and Dardanelles straits and our prediction of this impact should be consistent with geological records of relative sea-level (RSL) change in the vicinity of this area. As we have discussed, Lambeck et al. (2007) tuned their GIA simulation to be consistent with RSL records at four sites, the closest of which was 430–670 km from the two sills. Subsequent to their work, Vacchi et al. (2014) published RSL data from the northeast Aegean Sea and we have revised the ice history in the GIA simulation in order to fit the record from the site closest to the Bosphorus and Dardanelles straits, the island of Samothrace (Fig. 1). The sea level data were obtained from a site on the western tip of the island (Vacchi et al., 2014). Specifically, the ICE-5G ice history was modified by shallowing the eustatic curve by ∼2 m at 8 ka, and tapering this perturbation to 0 m by 4 ka, in order to provide a better fit to the local data set. This revised ice history, in combination with the 3-D Earth model described above, yields a good fit to the sea-level record at Samothrace (Fig. 4, red lines). As a point of comparison, we also show in Fig. 4 an RSL prediction at Samothrace based on the same ice history but using an Earth model that varies with depth alone (blue line). The model is identical to the spherically averaged rheology of the 3-D Earth model. The prediction of this model exhibits a much shallower sea-level trend since 10 ka, and significantly misfits the geological record of sea-level change. We performed a sequence of simulations based on a large suite of 1-D Earth models in which we varied the lithospheric thickness, upper mantle viscosity and lower mantle viscosity and found that none of the predictions using 1-D
S.L. Goldberg et al. / Earth and Planetary Science Letters 452 (2016) 178–184
181
Fig. 3. Lateral variations in mantle viscosity at a depth of 200 km below the western Black Sea and Aegean Sea incorporated into the GIA predictions described in the text. The symbol ν3D denotes the viscosity adopted in the calculation, and ν1D is the spherically averaged value of viscosity at this depth (5 × 1020 Pa s). The inset displays the variation in lithospheric thickness across the same region, including the location of plate boundaries (blue lines), incorporated into the lithospheric model. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
Fig. 4. Predictions of relative sea level at Samothrace (Fig. 1) for different GIA simulations. Red line: RSL prediction based on the ice history and 3-D Earth model discussed in the text, but no flooding of the Black Sea. Dashed, dotted red lines: Same as solid red line, except that flooding of the Black Sea that increased its bathymetry by an average of 50 m and 100 m, respectively, is included in the calculation. Blue line: Same as solid red line, except that lateral variations in Earth structure are ignored. The observations all involve the marine bivalve Cerastoderma glaucum and are taken from Vacchi et al. (2014). Two older observations reported in that study are not included. These are based on marine organics, represent terrestrial limits only and have potentially large age errors (A. Rovere, personal communication, 2016). The RMS misfit between the 3-D calculations and the observations is a factor of ∼3.2 smaller than the misfit associated with the 1-D calculation. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
Earth models provided a reasonable fit to the observations (see supplement, Fig. S1), suggesting that incorporating lateral variations in Earth structure is necessary to fit the observed data. To explore this issue further, Fig. 5 shows the difference in the RSL predictions at 10 ka over southern Europe and north Africa between the 3-D and 1-D Earth models. These differences are significant, with peak amplitudes of ∼10 m. The 3-D prediction is clearly influenced by large-scale variability in mantle viscosity and
lithospheric thickness; however, the impact of plate boundary effects is also noticeable in localized structure that extends across the region. For example, the southern boundary of the Eurasian plate is evident in a localized perturbation in sea level that extends west across Greece, and along the northern boundaries of the Aegean and Anatolian Plates, passing in the vicinity of the Dardanelles and Bosphorus straits. Since the boundary is modeled as a low viscosity zone within the lithosphere, the 3-D GIA calculation permits the crust on either side of the boundary to move relatively independently in response to the water loading that extends across the Mediterranean and Aegean. This permits greater subsidence of the Aegean crust (and thus sea level rise) in response to post-LGM meltwater loading, relative to the continuous lithosphere of the 1-D earth model prediction, which explains the results in Fig. 4. Two other sea level records exist in the Aegean Sea, further from the sills, on the islands of Limnos and Iskele (Vacchi et al., 2014). The fit of the 3-D model to these curves is less satisfactory than at Samothrace (see supplement, Figure S2). However, the incorporation of 3-D structure progressively improves the fit relative to the 1-D model as one considers sites closer to the sill. Given the short-wavelength variability evident in Fig. 5, we have chosen to tune the model to fit the site closest to the Black Sea sills. The solid red and blue lines in Fig. 4 and Fig. S2 were computed with the Black Sea masked so that no meltwater was allowed to enter it during the deglaciation phase. This raises the issue as to whether the flooding of the Black Sea would impact the sealevel trend predicted at Samothrace, and whether our tuning of the regional GIA model should account for this. To consider this, we ran simulations in which we flooded the Black Sea with varying volumes of water and tracked the resulting regional sea-level change over the next 10 kyr. As an example, Fig. 6 shows the total sea-level change 10 kyr after a flooding event that increased the bathymetry of Black Sea level by an average of 100 m. While the predicted impact of the flood on sea level is significant at the Bosphorus strait and sill, it is relatively minor at the Dardanelles sill (∼1 m) and at Samothrace (∼0.5 m). This is verified in Fig. 4, where the dashed and dotted red lines represent the predicted RSL curves at Samothrace for scenarios in which the original GIA calculation (solid red line) was augmented to include a Black Sea flood
182
S.L. Goldberg et al. / Earth and Planetary Science Letters 452 (2016) 178–184
Fig. 5. Difference in the predicted RSL at 10 ka (m) due to GIA based on calculations that adopt the 3-D and 1-D viscoelastic Earth models discussed in the text (i.e., the 3-D minus 1-D RSL predictions).
Fig. 6. Predicted total change in sea level (m) 10 kyr after a flooding event that raised the average bathymetry of the Black Sea by 100 m. The color on the main figure saturates at 12 m; the maximum value is ∼16 m in the center of the Black Sea. The prediction adopts the 3-D viscoelastic Earth model described in the text.
at 10 ka equivalent to an average bathymetry change of 50 m and 100 m, respectively. We conclude that our tuning of the GIA model to fit the RSL record at Samothrace is not sensitive to the timing or magnitude of the flooding event. However, in contrast to the predictions at Samothrace, the prediction in Fig. 6 indicates that the flooding event had a significant impact on the sea level change at the Bosphorus sill, and must therefore be taken into account when estimating the timing or magnitude of the flood event. To quantify this sensitivity, Fig. 7 shows predictions of GIA-induced RSL at Bosphorus based on the ice history discussed above and including various scenarios for the magnitude (Fig. 7A) and timing (Fig. 7B) of the Black Sea flood
event. The scenarios in Fig. 7A span the range of flood magnitudes proposed in the literature, and it is clear that these predictions are sensitive to this uncertainty. Indeed, relative to the no flood case, the 100 m event considered in Fig. 7A leads to a ∼10 m greater rise in the sea-level change predicted at the Bosphorus sill from 10 ka to present. This increased sea-level rise is largely due to crustal subsidence at the sill driven by the incursion into the Black Sea and associated water loading. Our complete set of calculations is used to construct Fig. 8, which plots RSL change at the location of the sill since the flood event versus the time of the event for several different flood magnitudes (in units of mean change in Black Sea depth). This figure
S.L. Goldberg et al. / Earth and Planetary Science Letters 452 (2016) 178–184
183
Fig. 7. Predictions of RSL change at the Bosphorus sill due to GIA and the Black Sea flood event. (A) Computed RSL change in the case of a flood event at 10 ka of varying magnitude (as labeled). (B) Computed RSL for a flood event of variable timing (as labeled) that raised the average depth of the Black Sea by 100 m. All calculations adopt the 3-D viscoelastic Earth model discussed in the text (see Fig. 3).
The results in Fig. 8 allow us to reassess previously published scenarios for the Black Sea flood event. For example, based on a synthesis of micropaleontological reconstructions, Soulet et al. (2011) estimated that the flood event occurred at 9 ka and, using GIA results in Lambeck et al. (2007), concluded that the sill has experienced 10 m of erosion during, and subsequent to, the flooding. Our results suggest that for a flood at 9 ka, erosion of 14–21 m occurred to bring the sill height to its present level of −33 m. If instead we adopt the 7.1 ka timing of Ryan et al. (1997), we estimate 24–30 m of erosion of the sill (i.e. if the sill had not eroded, it would be close to present sea level). In both cases, the range in the estimate reflects the uncertainty in the magnitude of the flood. 4. Discussion and conclusions
Fig. 8. Predicted RSL at the location of the Bosphorus sill at the time of imposed flooding for three different cases of flood magnitude (as labeled). Predictions are based on the ice history and 3-D Earth model described in the text.
allows for the estimation of one of the flood parameters (timing, magnitude, or sill depth) as a function of the other two. Consider, as an example, a flood event at 10 ka of magnitude 100 m. In this case, the RSL change since the flood is predicted to be −25 m. This prediction changes to −20.5 m and −17 m for flood events with the same timing and of magnitude 60 m and 30 m, respectively. Since the present day depth of the sill is ∼33 m, these scenarios imply that erosion of the sill over the last 10 kyr has lowered its elevation by 8 m, 12.5 m, and 16 m, for flood magnitudes of 100 m, 60 m, and 30 m, respectively. More generally, any scenario in Fig. 8 that yields a RSL above −33 m implies that the sill has eroded down to the observed level of −33 m since the time of flooding, while any scenario with a predicted RSL below −33 m implies that sedimentation has raised the elevation of the sill over the same time period. Thus, the results in Fig. 8 indicate that the sill must have experienced erosion if the flood event occurred at, or after, ∼11 ka. If instead the sill depth has not changed since the flood event, the event occurred between 11 and 11.6 ka.
We have predicted relative sea-level change at the location of the Bosphorus sill and surrounding regions since LGM using a gravitationally self-consistent ice age sea-level theory and a GIA model tuned to fit a RSL record in the northern Aegean. Our results indicate that the numerical sea-level predictions are sensitive to the incorporation of 3-D viscoelastic Earth structure, including lateral variations in mantle viscosity, lithospheric thickness and plate boundaries, as well as to the magnitude and timing of the flooding of the Black Sea. Neither of these sensitivities, which can perturb the sea-level predictions by as much as ∼10 m, was considered in the only previous GIA modeling of sea-level change in the region. We note that our interpretation of the RSL predictions has assumed tectonic stability over the time interval of study. While the area surrounding the Bosphorus is tectonically active, the magnitude of vertical tectonic motion that may have occurred over the past 10 ka is thought to be small. In the western Marmara Sea, which is a more tectonically active area, late Quaternary uplift rates have been calculated at ∼0.3 mm/yr, generally decreasing eastwards towards the Bosphorus (Yaltirak et al., 2002). Close to the Bosphorus itself, late Miocene limestones have been uplifted to 200 m, which indicates an average uplift rate of ∼0.02 mm/yr (Gökasan et al., 1997). This rate suggests an uplift of ∼0.2 m since 10 ka assuming no recent acceleration of the uplift, which is well within the uncertainties of our predictions. We have summarized our predictions in a plot showing RSL change at the location of the Bosphorus sill as a function of the
184
S.L. Goldberg et al. / Earth and Planetary Science Letters 452 (2016) 178–184
timing of the Black Sea flood for a suite of different flood magnitudes ranging from 30–100 m. As we noted in the Introduction, a wide range of estimates exist for the timing of the Black Sea flood and the issue has been the subject of significant, ongoing controversy – a debate that no doubt reflects the possible connection between the flood event and the Biblical flood legend. Our goal in the present study has been to provide a state-of-the-art geophysical analysis of this event that explicitly accounts for existing uncertainties in the nature of the sill and lake elevation at the time of the possibly cataclysmic event. Acknowledgements We thank two anonymous reviewers for constructive comments on a previous version of the manuscript. This work was supported by a HUCE summer research grant (to SLG), Harvard University, and the Canadian Institute for Advanced Research. Appendix A. Supplementary material Supplementary material related to this article can be found online at http://dx.doi.org/10.1016/j.epsl.2016.06.016. References Aksu, A.E., Hiscott, R.N., Yasar, D., Isler, F.I., Marsh, S., 2002. Seismic stratigraphy of Late Quaternary deposits from the southwestern Black Sea shelf: evidence for non-catastrophic variations in sea-level during the last ∼10000 years. Mar. Geol. 190, 61–94. Algan, O., Cagatay, N., Tchepalyga, A., Ongan, D., Eastoe, C., Gokasan, E., 2001. Stratigraphy of the sediment infill in Bosphorus Strait: water exchange between the Black and Mediterranean Seas during the last glacial Holocene. Geo Mar. Lett. 20, 209–218. Clark, P., Dyke, A., Shakun, J., Carlson, A., Clark, J., Wohlfarth, B., Mitrovica, J., Hostetler, S., McCabe, A.M., 2009. The Last Glacial Maximum. Science 325, 710–714. Conrad, C.P., Lithgow-Bertelloni, C., 2006. Influence of continental roots and asthenosphere on plate–mantle coupling. Geophys. Res. Lett. 33, L05312. Fichtner, A., Saygin, E., Tayaz, T., Cupilard, P., Capdeville, Y., Trampert, J., 2013. The deep structure of the North Anatolian Fault Zone. Earth Planet. Sci. Lett. 373, 109–117. Giosan, L., Filip, F., Constatinescu, S., 2009. Was the Black Sea catastrophically flooded in the early Holocene? Quat. Sci. Rev. 28, 1–6. Gökasan, E., Demirbag, E., Oktay, F., Ecevitoglu, B., Simsek, M., Yuce, H., 1997. On the origin of the Bosphorus. Mar. Geol. 140, 183–199. Hiscott, R.N., Aksu, A.E., 2002. Late Quaternary history of the Marmara Sea and Black Sea from high-resolution seismic and gravity-core studies. Mar. Geol. 190, 261–282. Kendall, R.A., Mitrovica, J.X., Milne, G.A., 2004. On post-glacial sea level – II. Numerical formulation and comparative results on spherically symmetric models. Geophys. J. Int. 161, 679–706.
Lambeck, K., Sivan, D., Purcell, A., 2007. Timing of the last Mediterranean Sea–Black Sea connection from isostatic models and regional sea-level data. In: YankoHombach, V.V., Gilbert, A.S., Panin, N., Dolukhanov, P.M. (Eds.), The Black Sea Flood Question: Change in Coastline, Climate and Human Settlement. Springer, Heidelberg, pp. 797–808. Lambeck, K., Smither, C., Johnston, P., 1998. Sea-level change, glacial rebound and mantle viscosity for northern Europe. Geophys. J. Int. 134, 102–144. Latychev, K., Mitrovica, J., Tromp, J., Tamisiea, M., Komatitsch, D., Christara, C., 2005. Glacial isostatic adjustment on 3-D Earth models: a finite-volume formulation. Geophys. J. Int. 161, 421–444. Lericolais, G., Guichard, F., Morigi, C., Minereau, A., Popescu, I., Radan, S., 2010. A post Younger Dryas Black Sea regression identified from sequence stratigraphy correlated to core analysis and dating. Quat. Int. 225 (2), 199–209. Major, C., Ryan, W., Lericolais, G., Hajdas, I., 2002. Constraints on Black Sea outflow to the Sea of Marmara during the last glacial–interglacial transition. Mar. Geol. 190, 19–34. Mitrovica, J., Forte, A., 2004. A new inference of mantle viscosity base upon joint inversion of convection and glacial isostatic adjustment data. Earth Planet. Sci. Lett. 225, 177–189. Nicholas, W.A., Chivas, A.R., Murray-Wallace, C.V., Fink, D., 2011. Prompt transgression and gradual salinisation of the Black Sea during the early Holocene constrained by amino acid racemization and radiocarbon dating. Quat. Sci. Rev. 30, 3769–3790. Peltier, W., 2004. Global glacial isostasy and the surface of the ice-age Earth: the ICE-5G (VM2) model and GRACE. Annu. Rev. Earth Planet. Sci. 32, 111–149. Ritsema, J., Deuss, A., van Heijst, H., Woodhouse, J., 2011. S40RTS: a degree-40 shear-velocity model for the mantle from new Rayleigh wave dispersion, teleseismic traveltime and normal-mode splitting function measurements. Geophys. J. Int. 184 (3), 1223–1236. Ryan, W.B.F., 2007. Status of the Black Sea flood hypothesis. In: Yanko-Hombach, V., Gilbert, A.S., Panin, N., Dolukhanov, P.M. (Eds.), The Black Sea Flood Question: Changes in Coastline, Climate, and Human Settlement. Springer, Dordrecht, pp. 3–88. Ryan, W.B.F., Pitman III, W.C., Major, C.O., Shimkus, K., Moskalenko, V., Jones, G.A., Dimitrov, P., Görür, N., Sakinç, M., Yüce, H., 1997. An abrupt drowning of the Black Sea shelf. Mar. Geol. 138, 119–126. Ryan, W., Pitman, W., 1998. Noah’s Flood: The New Scientific Discoveries About Events That Changed History. Simon and Schuster, New York. Ryan, W., Major, C., Lericolais, Gl., Goldstein, S., 2003. Catastrophic flooding of the Black Sea. Annu. Rev. Earth Planet. Sci. 31, 525–554. Soulet, G., Ménot, G., Lericolais, G., Bard, E., 2011. A revised calendar age for the last reconnection of the Black Sea to the global ocean. Quat. Sci. Rev. 30, 1019–1026. Spada, G., Alfonsi, L., Boschi, E., 1999. Chandler wobble excitation by catastrophic flooding of the Black Sea. Ann. Geofis. 42, 749–754. Vacchi, M., Rovere, A., Chatzipetros, A., Zouros, N., Firpo, M., 2014. An updated database of Holocene relative sea level changes in NE Aegean Sea. Quat. Int. 328–329, 301–310. Yaltirak, C., Sakinc, M., Aksu, A., Hiscott, R., Galleb, B., Ulgen, U., 2002. Late Pleistocene uplift history along the southwestern Marmara Sea determined from raised coastal deposits and global sea-level variations. Mar. Geol. 190, 283–305. Yanko-Hombach, V., Mudie, P., Kadurin, S., Larchenkov, E., 2014. Holocene marine transgression in the Black Sea: new evidence from the northwestern Black Sea shelf. Quat. Int. 345, 100–118.