Improved modelling of the Quaternary evolution of the Gulf of Corinth, incorporating erosion and sedimentation coupled by lower-crustal flow

Tectonophysics 440 (2007) 67 – 84 www.elsevier.com/locate/tecto

Improved modelling of the Quaternary evolution of the Gulf of Corinth, incorporating erosion and sedimentation coupled by lower-crustal flow Rob Westaway ⁎ Faculty of Mathematics and Computing, The Open University, Eldon House, Gosforth, Newcastle-upon-Tyne NE3 3PW, UK Received 11 January 2006 Available online 16 February 2007

Abstract Numerical modelling, incorporating coupling between surface processes and induced flow in the lower continental crust, is used to address the Quaternary evolution of the Gulf of Corinth region in central Greece. The post-Early Pleistocene marine depocentre beneath this Gulf overlies the northern margin of an older (Early Pleistocene and earlier) lacustrine basin, the Proto Gulf of Corinth Basin or PGCB. In the late Early Pleistocene, relief in this region was minimal but, subsequently, dramatic relief has developed, involving the creation of ∼ 900 m of bathymetry within the Gulf and the uplift by many hundreds of metres of the part of the PGCB, south of the modern Gulf, which forms the Gulf's main sediment supply. It is assumed that, as a result of climate change around 0.9 Ma, erosion of this sediment source region and re-deposition of this material within the Gulf began, both processes occurring at spatial average rates of ∼ 0.2 mm a− 1. Modelling of the resulting isostatic response indicates that the local effective viscosity of the lower crust is ∼4 × 1019 Pa s, indicating a Moho temperature of ∼ 560 °C. It predicts that the ∼ 10 mm a− 1 of extension across this ∼ 70 km wide model region, at an extensional strain rate of ∼ 0.15 Ma− 1, is partitioned with ∼ 3 mm a− 1 across the sediment source, ∼2 mm a− 1 across the depocentre, and ∼5 mm a− 1 across the ‘hinge zone’ in between, the latter value being an estimate of the extension rate on normal faults forming the major topographic escarpment at the southern margin of the Gulf. This modelling confirms the view, suggested previously, that coupling between this depocentre and sediment source by lower-crustal flow can explain the dramatic development in local relief since the late Early Pleistocene. The effective viscosity of the lower crust in this region is not particularly low; the strong coupling interpreted between the sediment source and depocentre results instead from their close proximity. In detail, the effective viscosity of the lower crust is expected to decrease northward across this model region, due to the northward increase in exposure of the base of the continental lithosphere to the asthenosphere; in the south the two are separated by the subducting Hellenic slab. The isostatic consequences of such a lateral variation in viscosity provide a natural explanation for why, since ∼ 0.9 Ma, the modern Gulf has developed asymmetrically over the northern part of the PGCB, leaving the rest of the PGCB to act as its sediment source. © 2007 Elsevier B.V. All rights reserved. Keywords: Greece; Corinth; Quaternary; Extension; Uplift; Subsidence; Lower-crustal flow

⁎ Also at: School of Civil Engineering and Geosciences, Newcastle University, Newcastle-upon-Tyne NE1 7RU, UK. E-mail address: [email protected]. 0040-1951/$ - see front matter © 2007 Elsevier B.V. All rights reserved. doi:10.1016/j.tecto.2007.02.002

68

R. Westaway / Tectonophysics 440 (2007) 67–84

1. Introduction The Gulf of Corinth in central Greece (Fig. 1) is one of the most rapidly extending active normal fault zones in the world and thus provides a superb environment for the study of continental extension. For most of its length, the south coast of this Gulf is bounded by northdipping normal fault segments (Fig. 2). Additional north-dipping active normal faults are present beneath the Gulf (e.g., Clement et al., 2004), and others are located in the northern part of the Peloponnese region to the south of it (e.g., Doutsos and Piper, 1990; Stewart and Vita-Finzi, 1996; McNeill and Collier, 2004). This Gulf is bounded to the south by an uplifted and partially eroded terrestrial (lacustrine) sedimentary basin (Fig. 2), known as the Proto-Gulf of Corinth Basin or PGCB (e.g., Ori, 1989). The southern part of the modern Gulf is underlain by the northern part of this basin. The youngest part of this lacustrine sequence has been dated to the late Early Pleistocene from mammalian biostratigraphy (see Ori, 1989, or Westaway, 1996, for a summary of evidence). Since deposition of this sequence ceased, the part of it beneath the modern Gulf has subsided dramatically, being now found beneath many hundreds of metres of Middle–Late Pleistocene sediment and up to ∼ 900 m of water. In contrast, much of the PGCB south of the modern Gulf has uplifted many hundreds of metres above sea-level. The Middle–Late Pleistocene succession that overlies the PGCB sediment beneath the Gulf largely consists of material eroded and reworked from the subaerial part of the PGCB to the south of the Gulf. It is now apparent that the present phase of extension in the Aegean region, associated with the modern geometry of the North Anatolian Fault Zone (Fig. 1), began in the Early Pliocene (e.g., Westaway, 2006b; Westaway et al., 2006a). There is thus no direct reason in terms of plate tectonics why the stable depositional environment represented by the PGCB should have been so dramatically disrupted in the past million years (cf. Armijo et al., 1996). The dramatic uplift of the northern Peloponnese region since the Early Pleistocene has thus been seen as anomalous. Many studies have indeed proposed ad hoc explanations for this “anomalous” footwall uplift, including crustal thickening caused by local crustal shortening (Mariolakos and Stiros, 1987), magmatism (Collier, 1990), and the isostatic response to the subduction of sediment within the Hellenic subduction zone (e.g., Le Pichon and Angelier, 1981). Armijo et al. (1996) suggested that this effect might be due to crustal thickening caused by a component of westward inflow of plastic lower crust from

beneath the subsiding central Aegean region. Westaway (1996, 1998) suggested instead that the main effect in this locality has been southward lower-crustal flow from beneath the modern Gulf interior to beneath the PGCB, driven by pressure variations at the base of the brittle layer caused by the sediment loading in the Gulf and the unloading of the PGCB as a result of its erosion. Westaway (2002a) followed this work up with quantitative numerical modelling of this combination of coupled processes, showing them to be a consequence of increased rates of surface processes following climate change at ∼ 900 ka. Marine oxygen isotope stage (MIS) 22, at ∼ 900 ka, is indeed well known to mark the first really large northern hemisphere glaciation, at the onset of predominant ∼ 100 ka (rather than ∼ 40 ka) Milankovitch forcing of climate, comparable in scale to the later large Middle and Late Pleistocene glaciations (e.g., Mudelsee and Schulz, 1997). By analogy with the last glacial cycle, the climate in the eastern Mediterranean region during MIS 22 can be anticipated to have been periglacial, with reduced vegetation cover resulting in increased rates of erosion (e.g., Westaway et al., 2006b). The Westaway (2002a) study was the first to apply this class of coupled model. It has led to similar studies of other regions, such as the analyses by Westaway et al. (2004, 2006b) of eroding Late Cenozoic terrestrial basins in western Turkey (Fig. 1) and the investigation by Morley and Westaway (2006) of some of the ‘superdeep’ Cenozoic basins of southeast Asia, whose post-rift evolution has also evidently involved coupling between surface processes and the induced flow in the lower continental crust. As a result of experience so gained, the numerical technique has been modified slightly (see Westaway et al., 2006b, for details). In addition, since the Westaway (2002a) study was prepared, significant new evidence pertaining to the Quaternary evolution of the Gulf of Corinth has developed. It is thus now appropriate to address this study region again. 2. Observational evidence The evidence available to constrain numerical solutions for the evolution of the Gulf of Corinth was discussed at length by Westaway (2002a). The principal requirement of any such modelling solution is to explain in an internally consistent manner, first, the increase in bathymetry of the Gulf from effectively zero in the late Early Pleistocene to the present maximum of ∼900 m, and, second, the uplift south of the Gulf that is revealed by the famous marine terrace staircase there (e.g., Keraudren and Sorel, 1987; Fig. 3). Paired instances such as this, of adjacent uplifting and subsiding areas

R. Westaway / Tectonophysics 440 (2007) 67–84

69

Fig. 1. Map of the Aegean region, showing active faulting and related sedimentation, adapted from Westaway (2006a, Fig. 1). The Gulf of Corinth is the sea area in the west-central part of the figure, between Patras and Athens, separating the land areas of Central Greece (Sterea Ellas) to the north and the Peloponnese (Peloponnisos) to the south. The island of Egina, a site of Quaternary hydrothermal activity associated with the subduction-related volcanism, is ∼ 15 km NE of Methana. The right-lateral North Anatolian fault zone enters the study region from the northeast, its strands terminating against northeast-dipping normal fault zones which bound the northeast coasts of Evvia and adjacent islands. Only Quaternary subduction-related volcanism is shown.

(the latter showing progressive increases in bathymetry, indicating that an excess of accommodation space over sediment supply is being dynamically created), are quite numerous, both in the eastern Mediterranean region and elsewhere worldwide; coupling between these areas, by sediment transport in one direction and induced lowercrustal flow in the opposite direction, is emerging as the general explanation for this effect (e.g., Westaway, 2004). The Gulf of Corinth can thus be regarded as a particularly dramatic example of this effect.

As in the Westaway (2002a) study, the section inland of Xylokastron (Figs. 2, 3) will be taken as representative of the uplift history for the area south of the central part of the Gulf. The oldest marine terrace in this area, which is thought to date from MIS 25 (see discussion of terrace ages by Westaway, 2002a), is found in this area at ∼ 700 m above sea-level, around Souli, ∼ 10 km south of the Gulf, indicating a time-averaged local uplift rate of ∼ 0.8 mm a− 1. The younger terraces, which crop out ∼ 1–2 km inland of Xylokastron, can be reliably

70

R. Westaway / Tectonophysics 440 (2007) 67–84

Fig. 2. Map of the Gulf of Corinth, adapted from Westaway (1996, Fig. 1). Map shows outcrop of PGCB lacustrine sediment (dot ornament) and Middle and Late Pleistocene marine terrace deposits (triangle ornament) south of the Gulf. Shading indicates offshore areas, with darker shading marking the ∼ 900 m deep abyssal plain beneath the central Gulf of Corinth and its smaller and shallower counterpart in the Alkyonides Gulf. Dashed line indicates the limit of the present-day drainage catchment south of the Gulf, which roughly matches the extent of the PGCB.

assigned to more recent MIS stages on the basis of the available absolute dating (also reviewed by Westaway, 2002a) and indicate uplift rates of ∼1.5 mm a− 1. However, as discussed by Armijo et al. (1996) and Westaway (2002a), these rates most likely exceed the rate at which localities farther south (e.g., around Souli) have uplifted on the same time-scale, due to the effect of back-tilting in the footwall of the major normal fault that bounds the coastline at Xylokastron (Fig. 3). The principal new data to emerge since the Westaway (2002a) study was prepared are as follows. First, it is now clear from several GPS publications (e.g., McClusky et al., 2000) that the extension rate across the Gulf of Corinth is ∼ 10 mm a− 1. The analysis for the Westaway (2002a) study pre-dated the McClusky et al. (2000) publication. On the basis of the geodetic information then available, which was much less clear, such high extension rates were considered, but the preferred extension rate adopted for the modelling was less; the modelling assumed an extensional strain rate in the crust of 0.06 Ma− 1, corresponding to an extension rate of ∼4 mm a− 1 across the 60 km wide extending

model region. However, such a low extension rate is no longer tenable. Second, studies of seismic wave propagation have improved constraint on the thicknesses of the crust and mantle lithosphere beneath the study region. Wide-angle seismic reflection results by Clement et al. (2004) indicate that the crustal thickness just north of the eastern Gulf is ∼ 32 km; gravity modelling by Tiberi et al. (2000, 2001) suggests that it is ∼ 31 km beneath the eastern Gulf but as low as ∼ 27 km beneath the PGCB farther south. I will assume that the crustal thickness at 0.9 Ma was 30 km throughout the model region. Tomographic studies by Tiberi et al. (2000) likewise indicate that the lithosphere beneath the Gulf of Corinth can be no more than ∼ 75 km thick, as the upper surface of the subducting Hellenic slab is imaged around that depth. The continental lithosphere in this area is thus somewhat thinner than the ∼ 100 km thickness assumed by Westaway (2002a). Third, the limited existing data on heat flow around the Gulf of Corinth (e.g., Fytikas and Kolios, 1979) has been supplemented by a new deep geothermal borehole

R. Westaway / Tectonophysics 440 (2007) 67–84

71

Fig. 3. Map of the south-eastern part of the Gulf of Corinth, adapted from Westaway (1996, Fig. 2) and based originally on Keraudren and Sorel (1987, Fig. 2) showing marine terrace localities (A to F) investigated by Westaway (1996). Terraces are numbered using the Keraudren and Sorel (1987) MIS assignments, although some of these ages are no longer preferred. Ornament indicates terrace correlations supported by Westaway (2002a).

at Egion (e.g., Rettenmaier et al., 2005; Fig. 2). Although results for it have yet to be published in detail, a temperature of 32 °C has been reported at ∼ 1000 m depth. Taking the mean surface temperature in this area as ∼ 12 °C, and the thermal conductivity, k, of the rocks that have been drilled as 2.5 W m− 1 °C− 1, heat flow of ∼ 50 mW m− 2 is indicated. Fytikas and Kolios (1979) inferred a gradual eastward increase in heat flow along the Gulf to ∼ 60 mW m− 2 in the Xylokastron area, followed by a more abrupt increase farther east. In their scheme, the surface heat flow increases to ∼80 mW m− 2 at Corinth town, then ∼ 90 mW m− 2 at Loutraki (which is of course a thermal resort with hot springs), before peaking well above 100 mW m− 2 in the Sousaki/ Egina/Methana area farther east (Figs. 1, 2), where hot springs, solfatara, and Quaternary volcanic activity are apparent. However, there was (and still is, as far as I am aware) no data to constrain the detail of this eastward increase, which could instead be more gradual, with values of maybe ∼ 70–80 mW m− 2 in the Xylokastron area. Fourth, seismic reflection profiling by Clement et al. (2004) has revealed significant active normal fault

escarpments beneath the Gulf of Corinth. It is thus evident that extension in this region is partitioned across faults in this location, as well as other active normal faults that were already known, bounding the coastlines and located to the south of the Gulf. This distribution of extension thus needs to be taken into account in future numerical modelling. Finally, considerable progress has been made relating deformation of the brittle upper crust in the Aegean region to changes in the geometry of subduction beneath it. Le Pourhiet et al. (2003) have suggested that the position of the Gulf of Corinth (as imaged by seismic tomography such as that by Tiberi et al., 2000), more or less directly above the point where the Hellenic slab starts to plunge steeply, means that this slab has significantly influenced the local stress field and thus controlled the geometry of normal faulting. The gradual transition northeastward from continental lithosphere that is shielded from the asthenosphere by the subducting slab to lithosphere that is exposed to the asthenosphere (and is indeed being heated by subduction-related volcanism in the Methana area) is presumably the reason why the surface heat flow varies so dramatically in the

72

R. Westaway / Tectonophysics 440 (2007) 67–84

vicinity of the Gulf of Corinth, as was noted above. Using a complex rheological model (few details of which were explained) Le Pourhiet et al. (2003) suggested that the uplift south of the Gulf is the isostatic response to local heating of the crust caused by southward retreat of the point at which the subducting slab starts to dip steeply. In other words, the crust south of the Gulf has begun to deform anomalously rapidly in the recent geological past because it has weakened due to the increased heat flow that has resulted because the continental lithosphere in this area has recently become exposed directly to the asthenosphere. However, a major problem with this interpretation is provided by the GPS dataset (e.g., McClusky et al., 2000), which indicates that the Gulf of Corinth is moving SSW at N 30 mm a− 1 relative to Eurasia, and thus probably moving even faster in this direction relative to the slab. It follows that the point where the slab begins to dip steeply is moving northward (or northnortheastward) relative to the Gulf of Corinth, as others (e.g., Laigle et al., 2004) have previously inferred, not southward as is required to drive the Le Pourhiet et al. (2003) numerical model. However, the effect of lateral variations in the thermal state of the crust, caused by the geometry of the subducting slab, on the numerical solutions presented in this study will later be discussed. 3. Modelling technique and results 3.1. Technique The modelling technique to be used is that developed by Westaway (2002a). It was fully defined in that study, including the underlying physics and the associated algebra, which included an analytic non-steady-state solution of the advective diffusion equation for heat transport within the crust, so these details are not repeated here. Minor modifications have since been made to the computer program, as documented by Westaway et al. (2006b). The basis of this technique, as applied to the present study region, is straightforward: erosion from the PGCB, south of the modern Gulf, unloads the underlying crust, reducing the local pressure at the base of the upper-crustal brittle layer. Likewise, loading of the crust beneath the Gulf, by sediment and seawater, increases the local pressure at the base of the brittle upper crust. The resulting lateral variation in pressure at the base of the brittle layer, Δp, acts to drive mobile lower-crustal material from beneath the depocentre to beneath the eroding sediment source, further affecting the overall isostatic balance. The model is treated as a

closed system for movement of crustal material: there is assumed to be no movement of any crustal mass into or out of the ends of the model or perpendicular to the plane of the model cross-section (Fig. 4). The only net change in mass is caused by any influx of seawater into the model should the Earth's surface in any part of it subside below sea-level. The start of each model run is marked by a step change in the rate of surface processes in each part of the model region. To facilitate comparison, this aspect has been kept the same as in the Westaway (2002a) study: before 0.9 Ma the whole model region, representing the PGCB, is assumed to have been a depocentre with a spatial average sedimentation rate of 0.1 mm a− 1, and that this sedimentation rate had persisted for long enough beforehand that this depocentre was in a thermal steady state. After 0.9 Ma, in response to effects of climate change, the model PGCB south of the modern Gulf began eroding at a spatial average rate of 0.2 mm a− 1, and the modern Gulf became a depocentre with a spatial average sedimentation rate also of 0.2 mm a− 1. For each 5 ka time step following these imposed changes, the computer program calculates the change in the pressure at the base of the brittle layer beneath each part of the model, taking account of the loss or gain of mass by surface processes, the vertical advection of the temperature threshold marking the base of the brittle layer, in response to these processes (given the nonsteady thermal state of the crust that results), and the thinning of the brittle upper crust due to the local extension. The program then uses the lateral pressure gradient that exists in each time step to calculate the resulting flow of lower crust between parts of the model.

Fig. 4. Summary of the model. The sediment source, hinge zone, and depocentre represent, respectively, the eroding part of the PGCB, the normal fault zone bounding the south coast of the Gulf of Corinth, and the depocentre beneath the Gulf. See text for discussion.

R. Westaway / Tectonophysics 440 (2007) 67–84

If any part of the model land surface falls below sealevel, it is assumed to become loaded by an appropriate height of water column, which thus also contributes to Δp. The mantle lithosphere is assumed to be extending along with the crust, and to thus be thinning. The isostatic consequences of this thinning are incorporated into the calculations of vertical crustal motions. However, since the initial width of the extending model region (set to 60 km, as for Westaway, 2002a) is small compared with the likely flexural half-wavelength for the mantle lithosphere, no warping of the mantle lithosphere is permitted; this layer is assumed instead to experience uniform vertical adjustments in response to the spatial average load applied to it. In order to distinguish effects of extension on normal faults within the eroding part of the PGCB and beneath the floor of the modern Gulf from those of other normal faults that delineate the model ‘hinge zone’ (Fig. 4) where the land surface drops between these model regions, the computer program allows separate extensional strain rates to be specified: ET for the model crust as a whole, Ee for the sediment source region and Es for the depocentre (Fig. 4). As in the Westaway (2002a) study, the program also allows the extensional strain rate in the mantle lithosphere, Em, to be specified independently of that in the crust. This facility is provided to enable investigation of the consequences of ‘depth-dependent stretching’, a concept that features in many analyses of extensional sedimentary basins by others. However, this facility will not be used in this study; all solutions will assume that Em = ET. Calculations of rates of flow in the lower crust assume that this layer is isoviscous, with effective viscosity ηe. The interpretation of ηe in terms of the real profile of viscosity decreasing with depth due to the increase in temperature follows Westaway (1998). This calibration is appropriate for a temperature-dependent linear viscous rheology for a material with a uniform geothermal gradient in the lower crust, with a molar activation energy G = 200 kJ mol− 1, and with a temperature of 350 °C marking the brittle–ductile transition at the top of this viscous layer. In other words, it is assumed that viscosity η varies with absolute temperature T as η(T) = A × exp (−G / R T). The scale factor A is taken as 200 kPa s− 1, after Westaway (1998), and the molar gas constant R, as 8.3143 J mol− 1 °C− 1. As noted by Westaway (1998), the viscosity of the deepest crust, just above the Moho, ηm, can be estimated using this approach as c. ηe / 60; the associated Moho temperature, Tm, is then calculated using the above viscosity calibration. To illustrate the effect of

73

different assumed values of ηe on the isostatic response, each of the solutions presented in this study will be run for 6 different values of ηe, listed along with the corresponding Moho temperature calculated using the above calibration: M1, ηe = 8 × 1019 Pa s; Tm = 542 °C; M2, ηe =4×1019 Pa s; Tm =561 °C; M3, ηe =2×1019 Pa s; Tm =582 °C; M4, ηe =1× 1019 Pa s; Tm = 603 °C; M5, ηe =5×1018 Pa s; Tm =626 °C; and M6, ηe =2.5×1018 Pa s; Tm =650 °C. Graphs of viscosity versus depth and of steady-state flow velocity versus depth for this form of viscosity calibration have been published previously, by Westaway (1998, 2002b, in press-b). The graphs in Fig. 7b of Westaway (2002b) and Fig. 1b of Westaway (2006c) illustrate conditions similar to those deduced in the present study for the lower crust beneath the Gulf of Corinth, and so are not repeated here. As previous studies (e.g., Westaway, 2002a; Morley and Westaway, 2006) have shown, this technique allows the display of variations over time of many different parameters derived from the model. For the purpose of this study, the output will be restricted to graphs of predicted uplift and uplift rate of the sediment source and of bathymetry within the Gulf, since each of these can be compared with observational evidence. By ‘uplift’ and ‘uplift rate’ I mean the increase in altitude (and its rate) of markers that are not themselves eroding, such as the surfaces of marine terraces, not the envelope of successive typical altitudes of the eroding land surface. In addition, graphs of the parameter Δp / ηe will be shown. As already noted, the technique requires Δp N 0 to drive lower-crustal flow in the sense described above. The purpose of displaying this parameter is twofold. First, critics (Allen et al., 2004) have asserted (without providing any evidence) that the technique used in this study is wrong in principle, because it ‘ought’ in their view to predict lateral pressure gradients in the opposite sense. Explicit display of calculated values of Δp / ηe demonstrates, on the contrary, that this is not so. Second, it has been noted previously (e.g., Westaway et al., 2006b), following other modelling of coupled systems of this type, that adjustment of Δp / ηe back towards a new thermal and mass-flux steady state consistent with the model rates of erosion and sedimentation typically occurs in a manner that is independent of ηe (i.e., in the limit of long time-scales after the change in rates of surface processes, Δp ∝ ηe so Δp / ηe is a constant). This test can be applied in the present modelling to see how close the 0.9 million years that have elapsed since the rates of model surface processes have increased comes to re-establishing steadystate conditions.

74

R. Westaway / Tectonophysics 440 (2007) 67–84

3.2. Results Fig. 5 shows the first set of model results. To derive it, the initial (i.e., at 0.9 Ma) widths of the model sediment source and depocentre, Le and Ls have each been set to 30 km and Lh, the nominal width of the hinge zone in between, has been set to 4 km (all these values, like those of all other model parameters that are not discussed, are the same as were used by Westaway, 2002a). ET has been set to 0.2 Ma− 1, giving 12 mm a− 1 of extension across the 60 km width of model crust that is extending. The initial thicknesses of the crust and mantle lithosphere, Hco and Hmo, have been set to 30 km and 60 km. The initial altitude of the land surface

is set to sea-level, thus as soon as there is any net subsidence of the sediment surface in the depocentre, it will become seawater-loaded and this water load will add to the calculated subsidence. Unlike in the Westaway (2002a) study, no attempt is made to consider complications due to the fluctuating water level in the Gulf in response to glacio-eustatic sea-level variations (discussion of this point contributed significantly to the complexity of that explanation without really adding anything of substance). Having specified these values, along with the others mentioned earlier, the only freedom remaining is to vary Ee and Es in order to match observations as well as possible. Such matching needs to take account of the

Fig. 5. Set 1 of model solutions. (a) Predicted variations over time in the spatial average depth of the Gulf of Corinth. Solid symbol indicates the present-day maximum depth (∼900 m). The present-day spatial average depth across the modelled cross-section is probably 600–700 m. (b) Predicted variations over time in the spatial average uplift south of the Gulf of Corinth. Solid symbols indicate uplift estimates from reliably dated marine terraces. Open symbols indicate alternative age assignments or amounts of uplift for terraces whose ages are not reliably constrained. Terrace age and height data are from Westaway (2002a), which lists original sources of information. Terraces are labelled with MIS stage numbers at their times of formation, substages being not labelled to avoid clutter. The spatial average uplift across the modelled cross-section in the most recent part of the record is thought to be significantly less than that indicated by the youngest marine terraces, all of which are close to Xylokastron, immediately south of the south coast of the Gulf. (c) Predicted variations over time in the spatial average uplift rate south of the Gulf of Corinth. Solid symbols indicate time-averaged values between reliable data points in (b). The spatial average uplift rate across the modelled cross-section in the most recent part of the record is thought to be significantly less than that indicated by the youngest marine terraces, all of which are close to Xylokastron, immediately south of the south coast of the Gulf. (d) Predicted variation in the lateral pressure variation Δp at the base of the brittle upper crust between the Gulf of Corinth depocentre and the eroding sediment source farther south, normalised for each solution by ηe, the effective viscosity of the lower crust. Positive values of Δp, as shown, act to drive lower crust from beneath the depocentre to beneath the sediment source. See text for discussion.

R. Westaway / Tectonophysics 440 (2007) 67–84

fact that the model works with spatially-averaged parameter values that are uniform in cross-section (Fig. 4), whereas along any cross-section across the real Gulf observables such as bathymetry and uplift rates vary laterally. For instance, the maximum bathymetry in the Gulf ‘abyssal plain’ is ∼ 900 m, but this only occurs in a small area just north of the south coast. North of it, the Gulf is shallower, the bathymetry tapering gently northward towards the north coast. I thus estimate crudely that the mean bathymetry for a section through the central Gulf is maybe two thirds or three quarters of the maximum, or ∼ 600–700 m. Likewise, the southward back-tilting means that the uplift rates of the youngest marine terraces in the immediate vicinity of Xylokastron are maybe double the contemporaneous spatial-average rates for the PGCB south of the Gulf. Fig. 5 shows a set of six uplift modelling solutions, consistent with these constraints, for different values of ηe. The values of Ee and Es have been adjusted to best fit the observational constraints. The optimum solution, characterised by Ee = 0.010 Ma− 1 and Es = 0.065 Ma− 1, is roughly intermediate between predictions M2 and M3, thus indicating a Moho temperature of ∼ 570 °C (Table 1). This solution predicts extension at ∼ 3 mm a− 1 across the depocentre and ∼ 2 mm a− 1 across the sediment source, thus indicating extension at ∼ 7 mm a− 1, at the start of the model run, across the major normal fault zone bounding the south coast of the Gulf, which forms the ‘hinge zone’ in the model. The predicted mean present-day thicknesses for this solution are 25.1 km for the model crust and 50.1 km for the mantle lithosphere, making the present-day overall lithosphere thickness 75.2 km. Assuming an asthenosphere temperature of 1300 °C, the geothermal gradient

75

in the mantle lithosphere, co, can be determined as 14.6 °C km− 1 (i.e., 1300–570 °C/50.1 km), corresponding to a basal heat flow (calculated as co × k, with – as before – k = 2.5 W m− 1 °C− 1) of 37 mW m− 2. Extrapolating this geothermal gradient into the crust indicates that 204 °C of the temperature rise there is attributable to radioactive heating, rather than basal heat flow. Using relevant theory (from Lachenbruch, 1970; see also Westaway et al., 2006b), this radioactive contribution can be explained by appropriate combinations of Yo, the radioactive heat production in the uppermost crust, and D, the scale depth over which radioactive heating decreases to 1/e of its maximum value. Assuming D = 10 km, the radioactive contribution to the near-surface geothermal gradient, cr, would be 20.4 °C km − 1 , which when added to the basal contribution would make the near-surface geothermal gradient, cs, 35.0 °C km− 1. The surface heat flow qs can then be calculated as cs × k and is 87 mW m− 2. An equivalent calculation with D = 15 km would instead yield qs = 71 mW m− 2 (Table 1). The principal features of these model solutions are as follows. The model sediment source begins to uplift as soon as erosion of it begins. In contrast, the model depocentre remains initially at (or very close to) sealevel. After a time-lag (which depends on ηe, decreasing from 195 thousand years for solution M1 to 50 thousand years for solution M6) the model depocentre begins to subside below sea-level. The onset of water-loading at this time affects the overall isostatic response, causing an inflection in each of the graphs (most obvious in those of uplift rate and Δp / ηe in Fig. 5c and d). From this time onward, for several tens of thousands of years, the effect of water-loading on the model is small. This

Table 1 Predictions of extension rates and the thermal state of crust in the vicinity of the Gulf of Corinth D = 10 km

D = 15 km

Ee Es Veo Vso ηe Tm Hc Hm co cr qs cr qs Hmo. ET (km) (Ma− 1) (Ma− 1) (Ma− 1) (mm a− 1) (mm a− 1) (1018 Pa s) (°C) (km) (km) (°C km− 1) (°C km− 1) (mW m− 2) (°C km− 1) (mW m− 2) 60 55 50

0.200 0.150 0.150

0.100 0.100 0.105

0.065 0.065 0.070

3.0 3.0 3.2

2.0 2.0 2.1

30 40 40

570 25.1 50.1 14.6 560 26.2 48.2 15.4 560 26.2 43.7 16.9

20.4 15.7 11.7

87 78 71

13.6 10.5 7.8

71 65 62

Hmo is the mantle lithosphere thickness assumed at the start of each model run (i.e., at 0.9 Ma). ET is the assumed spatial average overall extensional strain rate across the model region. Ee and Es are the best-fitting spatial average extensional strain rates across the model sediment source and depocentre. Veo and Vso are the corresponding initial extension rates across the model sediment source and depocentre, calculated as Veo = Ee × Leo and Vso = Es × Lso, where Leo and Lso are the initial widths of the sediment source and depocentre (both 30 km). ηe and Tm are the predicted effective viscosity of the lower crust and Moho temperature, calculated as described in the text. Hc and Hm are the predicted present-day mean thickness of the crust and thickness of the mantle lithosphere. The method used for prediction of the present-day thermal state of the lithosphere is explained in the text: co is the basal geothermal gradient, cr the contribution to the near-surface geothermal gradient from radioactive heating in the crust, D the scale depth for this radioactive heating, and qs the surface heat flow. Calculations of cr and qs assume either D = 10 km or D = 15 km. Other parameters required to define the model are specified in the text and/or take the same values as were used by Westaway (2002a,b).

76

R. Westaway / Tectonophysics 440 (2007) 67–84

effect acts to perturb the model further away from steady state, whereas without it the model would begin much sooner to adjust back to a new steady state consistent with its new rates of erosion and sedimentation. At a time 110 thousand years after the start of solution M6, and later for the higher-viscosity solutions, by which point the model water depth has reached ∼ 60 m, the water-loading of the depocentre ‘wins’ this competition. From this time onward until the present day the progressive increase in water-loading drives lower crust out from beneath the Gulf at ever increasing rates, progressively creating the present dramatic relief across its southern margin. As Fig. 5d indicates, the model pressure gradient at the base of the brittle upper crust is southward, thus driving lower-crustal flow outward from beneath the depocentre to beneath the sediment source (cf. Allen et al., 2004). In the higher-viscosity solutions (M1 to M3) the mean geothermal gradient in the model upper crust remains at all times less beneath the depocentre than beneath the sediment source, thus making a net contribution of its own to the lateral pressure gradient, supplementing the larger contribution from the water load. In the lower-viscosity solutions (M4 to M6) the high mobility of the lower crust instead means that so

much of it flows out from beneath the Gulf during the early stages of each model run that in the later stages the brittle upper crust becomes thinner beneath the depocentre than beneath the sediment source. In these cases, the water load thus counteracts the tendency that would exist, if this load were removed, for this material to flow back to beneath the Gulf. The fact that in these solutions the water-loading and upper-crustal-thickness contributions to Δp are opposed is the main reason why, as ηe is reduced, the graphs of parameter variations in Fig. 5 converge (i.e., the predictions for lower values of ηe differ much less from each other than those for higher values of ηe, even though the ratio of ηe values between different solutions is constant). Further reductions in ηe would then make little difference to the overall isostatic response. Conversely, increasing ηe to values much higher than that for solution M1 would reduce the predicted amounts of lower-crustal flow, thus reducing the relief that is predicted to develop. In the limit of very high ηe the solution would tend towards that predicted for conventional Airy isostasy without lower-crustal flow, and would indeed coincide with syn-rift predictions of the McKenzie (1978) stretching model. However, for all values of ηe considered, the predicted values of Δp / ηe continue to diverge at the present day

Fig. 6. Set 2 of model solutions. Data sources and display format are the same as for Fig. 5. See text for discussion.

R. Westaway / Tectonophysics 440 (2007) 67–84

(Fig. 5d). This indicates that the model crust in this region continues to adjust away from a thermal and mass-flux steady state, as a result of the magnitude of the surface processes applied to it (principally, the magnitude of the water-load, as discussed above). Experience with applying this type of model on longer timescales (Westaway et al., 2006b; Morley and Westaway, 2006) indicates that it typically takes several million years, at least, for graphs of Δp / ηe for different values of ηe to converge, and even longer for such coupled systems to reach a new steady state consistent with their rates of surface processes. It is of course unlikely that this state will ever be reached in the Gulf of Corinth, because the sediment supply from the eroding part of the PGCB will be exhausted long beforehand. It is evident that even though ηe in this region is predicted to be not particularly low, the isostatic consequences of the induced lower-crustal flow are dramatic. As Westaway (2002a) pointed out, the very strong coupling, by lower-crustal flow, between this depocentre and sediment source results from their close proximity. The corresponding sediment transport in the opposite sense occurs via the many short, high-gradient rivers that drain from the PGCB into the southern margin of the Gulf, as illustrated in Fig. 2.

77

As well as fitting the observational evidence regarding vertical crustal motions, this solution complies with the geometrical constraint on present-day lithosphere thickness imposed by the geometry of the subducting slab (see above) and predicts a range of possible values for the surface heat flow that are in broad agreement with what is expected (see above). The principal problems with it are that the predicted mean crustal thickness is rather low and the predicted extension rate is rather high; having started 60 km wide at 0.9 Ma, the extending model region is predicted to now be 70.8 km wide, making its predicted present-day extension rate 14.2 mm a− 1. In an attempt to mitigate these problems, I now consider a second solution, in Fig. 6. Here, the extensional strain rate has been reduced to 0.15 Ma− 1, such that the predicted present-day width of the extending model region is 68.1 km and its present-day extension rate is 10.2 mm a− 1. Because the predicted thinning is slower than before, the initial thickness of the mantle lithosphere must be made less than before, no more than 55 km, for the present-day thickness of the lithosphere to fit the constraint from the geometry of subduction. The solution that best fits the observational evidence is now M2, indicating a predicted Moho temperature of

Fig. 7. Set 3 of model solutions. Data sources and display format are the same as for Fig. 5. See text for discussion.

78

R. Westaway / Tectonophysics 440 (2007) 67–84

∼ 560 °C. The predicted basal heat flow is thus slightly greater than before, 39 mW m− 2. This means that the radioactive contribution to surface heat flow is less than before, being 39 mW m− 2 for D = 10 km and 26 mW m− 2 for D = 15 km, making the predicted surface heat flow in the range 65–78 mW m− 2 (Table 1). In the third solution, in Fig. 7, I reduce the initial mantle lithosphere thickness further, to 50 km, keeping the extensional strain rate at 0.15 Ma− 1. The best-fitting solution is, once again, M2, again indicating a predicted Moho temperature of 560 °C. The predicted rates of extension within the depocentre and sediment source, rather than within the hinge zone, are now slightly greater than before (Table 1). The reduced mantle lithosphere thickness means that the present-day basal heat flow is greater than before, 42 mW m− 2. This means that the radioactive contribution to surface heat flow is less (Table 1), making the predicted surface heat flow in the range 62–71 mW m− 2. Other solutions are also possible. However, in general the trade-offs evident between Figs. 6 and 7 continue to apply: further reduction in the initial thickness of the mantle lithosphere, Hmo, requires further slight increases in rates of ‘distributed’ extension, but predicted best-fitting values of ηe and Moho temperature do not change significantly. These adjustments are thus accompanied by further predicted increases in basal heat flow and reductions in the predicted radioactive contribution. The limiting case occurs when Hmo is ∼ 40 km, for which the predicted radioactive heating would be zero. In this case the present-day lithosphere thickness would be ∼61 km and the surface heat flow would entirely comprise a basal contribution, which would be 53 mW m− 2. However, the surface heat flow in the part of the Gulf of Corinth that is being modelled is significantly greater than this threshold (see Fytikas and Kolios, 1979, also earlier discussion above), from which it follows that the lithospheric thickness in this region is somewhat in excess of 61 km. 4. Discussion The preferred estimates of ηe and Moho temperature obtained in this study, ∼ 4 × 1019 Pa s and ∼ 560 °C, are very similar to those derived by Westaway (2002a) for the same region using the same technique, which were ∼ 6 × 1019 Pa s and ∼550 °C. The principal difference between the two sets of results is that the new set is consistent with a much higher extension rate across the Gulf of Corinth, ∼10 mm a− 1 against ∼ 4 mm a− 1. As already noted, Westaway (2002a) was unable to find any

satisfactory model solution consistent with such a high extension rate. However, all the solutions considered in that study involved the same strain rates for distributed extension in the sediment source and in the depocentre (i.e., Ee = Es). Now that the computer program has been modified to allow solutions with these strain rates assigned different values, it has become clear that, if the overall extensional strain rate ET is set to relatively high values, the only viable solutions for this study region involve Ee somewhat greater than Es (Table 1). Increases in rates of vertical crustal motions, starting around the time of MIS 22, and attributed (like that in the present study) to non-steady-state isostatic consequences of climate change (such as increased rates of erosion during this and subsequent glacial maxima, when this region developed a periglacial climate), are now widely recognised in the eastern Mediterranean region (e.g., Bridgland et al., 2003; Demir et al., 2004; Westaway et al., 2004, 2006a,b; Westaway, 2006b; Bridgland et al., in press; Demir et al., in press; Westaway and Bridgland, in press). Currently, the ‘benchmark’, representing the most tightly-constrained analysis of this type, is the modelling undertaken to explain the post-Pliocene uplift of the Kula area of western Turkey (Fig. 1) (e.g., Westaway et al., 2006b). In this locality a former lake basin became disrupted by fluvial incision starting at ∼ 3 Ma; since this time, the Gediz River has locally incised by ∼ 400 m. This incision history is constrained by a remarkable staircase of river terraces, many of which have been dated by K–Ar and Ar–Ar dating of capping basalt flows (e.g., Westaway, 2004; Maddy et al., 2005; Westaway, 2006a) and can thus be assigned to individual MIS stages, even in the Early Pleistocene (providing a degree of age control that is virtually unmatched worldwide). Numerical modelling of this uplift history indicates that ηe is locally ∼ 2 × 1018 Pa s, making the Moho temperature ∼ 660 °C, indicating much hotter crust than in the vicinity of the Gulf of Corinth. Aside from the super-accurate age control, this part of western Turkey is much simpler in principle to model than Corinth, because there is no significant local extension, so complexities associated with constraining local strain rates can be excluded (the region is far enough from the Alaşehir Fault Zone, Fig. 1, for this structure not to have any measurable effect). In addition, the local surface heat flow (110–120 mW m− 2) has been accurately measured by applying modern hydrogeochemical techniques to hot springs (e.g., Ilkışık, 1995). As noted above, on the contrary, the limited evidence to constrain surface heat flow and radioactive heat production in and around the Gulf of Corinth is a significant factor limiting the degree of constraint now

R. Westaway / Tectonophysics 440 (2007) 67–84

possible in the modelling of this region. Furthermore, past changes to the geometry of subduction beneath western Turkey have a minimal effect on the thermal state of the lithosphere (Westaway, 2006a); being so thin (∼ 60 km thick), this layer retains little long-timescale ‘thermal memory’ of such events. In contrast, as already noted, the presence of the subducting Hellenic slab can be expected to introduce a lateral gradient in the thermal state of the crust across the Gulf of Corinth and its surroundings. One potential problem with the numerical technique used in the present study, noted previously (Morley and Westaway, 2006), is that even if ηe is initially the same everywhere within a model region, then as time

79

progresses, the lateral variations in crustal thickness that are induced by the model will result in lateral variations in Moho temperature and thus in ηe. However, the relatively limited lateral variations in crustal thickness predicted by the solutions in Figs 5, 6 and 7 will result in lateral variations in Moho temperature by no more than ∼ ± 10 °C relative to the mean value, with the deepest crust becoming systematically hotter beneath the sediment source and cooler beneath the depocentre. This induced effect will thus counteract any initial lateral variation in the thermal state of the crust due to the Hellenic slab; consequently, although the modelling program in its present form does not take into account either of these potential sources of systematic error, they

Fig. 8. Cartoon providing a schematic illustration to suggest why the modern Gulf of Corinth developed near the northern margin of the Proto Gulf of Corinth Basin at the end of the Early Pleistocene. (a) The stable, symmetrical, lacustrine depocentre envisaged for the PGCB in the late Early Pleistocene. Grey shading indicates the most recent sediments deposited. (b) Mass fluxes (in arbitrary units) following climate-induced reduction in lake level inferred during MIS 22. Black arrows denote directions and magnitudes of sediment fluxes; white arrows denote directions and magnitudes of fluxes of lower-crustal flow. S and E denote regions of sedimentation and erosion at this time. Dark grey shading indicates the sediment deposited at this stage. In response to the same magnitude of sediment fluxes across the basin, the fluxes of lower crust increase from south to north, consistent with the northward decrease in ηe and northward increase in Moho temperature that are inferred to have existed, given that the continental lithosphere in the north was exposed directly to the asthenosphere but that in the south was shielded from it by the subducting Hellenic slab. The relative magnitude of the expected lateral variation in fluxes of lower-crustal material has been exaggerated for emphasis. (c) Net mass fluxes into (+) or out of (−) each part of the diagram resulting from the processes in (b) and associated changes in down-to-the-north relief. These net mass fluxes determine the resulting overall change in the level of the land surface. It involves a decrease in surface altitude in the north of the PGCB, where outflow of lower crust exceeds the influx of sediment, and an increase in much of the south, where influx of sediment exceeds outflow of lower crust. As a result of these changes, the basin is no longer symmetrical; the lower land surface in the south is conducive to future sedimentation being concentrated there, whereas the increase in altitude in the south is conducive to future erosion from there, once an overall south–north topographic gradient has been established.

80

R. Westaway / Tectonophysics 440 (2007) 67–84

will to some extent cancel in this particular set of solutions. As previously noted, Le Pourhiet et al. (2003) tried to explain the observed dramatic asymmetry of the modern Gulf of Corinth. Their explanation is almost certainly wrong (see above), but they nonetheless raised valid questions, relating to the geometry of this Gulf. Why is it located across the northern margin of the PGCB? Why is it not centrally-located or at the southern margin of this older basin? Is this an effect of chance, or might it possibly be a consequence of the lateral variation in thermal state of the crust, caused by the Hellenic slab? Fig. 8 attempts to address this issue semi-quantitatively. Part (a) shows in cross-section the assumed form of the PGCB in the late Early Pleistocene, just before MIS 22. It is assumed to have been symmetrical, with sediment accumulating in a shallow lake, draping across normal fault footwalls and largely filling the topographic lows in hanging-walls. However, the asymmetry was partly broken by the northward decrease in ηe in the lower crust, expected from the presence of the Hellenic slab. Suppose that during MIS 22, climate change caused a reduction in lake level (whether as a result of sea-level fall due to glacio-eustasy, or reduced rainfall due to cooling of the seas and atmosphere at temperate latitudes, or whatever). During this cold stage there would have been little vegetation to stabilise slopes, so local footwall highs would be expected to begin to erode significantly, the eroded material accumulating in the adjacent hanging-wall lows (Fig. 8b). These local loading and unloading effects would be expected to induce local patterns of lower-crustal flow. However, the magnitude of the induced lower-crustal-flow response would be expected to vary laterally due to the lateral variation in ηe, being most vigorous in the north and least vigorous in the south (Fig. 8b). As Fig. 8c shows, the outflow of lower crust will thus be greatest from beneath the hanging-wall of the northernmost normal fault system within the PGCB. This is indeed the only place in the figure where absolute subsidence of any hanging-wall locality is predicted as a result of this process, and also the only place where additional down-to-the-north relief is created that adds to the existing relief (elsewhere, where relief in this sense is created, it partly cancels existing relief in the opposite sense, between hanging-wall and footwall localities). The effect of this combination of processes is thus to remove the initial symmetry by acting to shift the lowest part of the land surface to near the northern margin of the basin. A subsequent rise in the lake level would thus have the effect of making this northernmost part of the basin the principal sediment trap. Over time,

after cyclic repetition of this combination of processes, one would thus expect to see a gradual adjustment, whereby deposition became concentrated in the northern part of the former PGCB, with its southern part developing into an eroding landscape — as is observed today. However, it can be anticipated that, in detail, this transition involved a complex sequence of local adjustments, which took place during the initial tens or maybe hundreds of thousands of years following MIS 22. Evidence for this adjustment process may well thus exist in the basal part of the post-PGCB sedimentary succession within the modern Gulf, which has not yet been studied in detail. However, erosional unconformities are locally evident, on seismic reflection profiles (e.g., Clement et al., 2004), between the PGCB and postPGCB successions beneath the Gulf, consistent with the expected complex sequence of adjustments shown schematically in Fig. 8. This complex initial adjustment process may correspond to the initial part of each of the model solutions in Figs. 5, 6 and 7 where inflections occur (for instance in the predicted uplift rate), before the isostatic response within the modern Gulf depocentre became dominated by water-loading. As noted previously, the present modelling technique does not allow for lateral variations in rheology, so such effects cannot be investigated directly, although they can be investigated in a semi-quantitative way, as in Fig. 8. Likewise, the technique spatially averages both the extensional strain and the vertical crustal motions across each part of the model region. This makes it impossible to relate any model prediction to some forms of detailed local evidence, such as amounts of tilting of individual marine terraces and/or offsets of these terraces across individual normal faults. Other key observables available for the Gulf of Corinth are provided by gravity studies, geodesy and earthquake seismology. As already noted, the GPS evidence indicating ∼ 10 mm a− 1 of extension across the Gulf (e.g., McClusky et al., 2000) has been used as an observational constraint in the present modelling, and so cannot be used to test the model predictions. Many studies of microseismicity indicate that microearthquake focal depths persist to ∼ 15 km depth in this study region, thus marking the base of the brittle layer, at depth zb, that is taken in this study as the 350 °C isotherm. However, such focal depth determinations are typically uncertain by several kilometres, thus making completely pointless any attempt at using them to look for the small lateral variations in zb that are predicted by the modelling. For instance, for solution M2 in Fig 7 the present-day values of zb are predicted to be 14,149 m beneath the Gulf and 14,005 m beneath the eroding sediment source; similar small

R. Westaway / Tectonophysics 440 (2007) 67–84

differences in zb are predicted for the other solutions. There is no possibility of seismograph station coverage in this region being improved to the extent that would be necessary to make such an observational test possible. The most telling of these geophysical observables, for use in an observational test, is thus the gravity dataset. Tiberi et al. (2001) produced a Bouguer gravity map of the Gulf of Corinth and its surroundings, based on thorough compilation of the available data. To a first approximation, this map shows the Bouguer gravity anomaly decreasing westward from ∼ + 40 mgal at the eastern end of the Gulf to ∼− 100 mgal at its western end. Tiberi et al. (2001) explained this effect as a consequence of an east–west increase in crustal thickness. Superimposed on this variation is a much smaller-amplitude variation across the Gulf, the Bouguer gravity anomaly along any profile being typically ≤ 10 mgal less beneath the Gulf than beneath the northern Peloponnese farther south. This would correspond, using the parameters deduced from the Tiberi et al. (2001) analysis, to a typical crustal thickness ∼ 3 km greater beneath the Gulf than beneath the northern Peloponnese (∼ 31 km against ∼28 km for the eastern part of the Gulf. However, analysis of seismic wave propagation, also by Tiberi et al. (2001) indicated that the crust beneath parts of the northern Peloponnese may be as much as 2 km thicker than was predicted by their own gravity modelling. Tiberi et al. (2001) also estimated that the uncertainty in each of their gravity data were ∼ ± 4 mgal. Taking all these uncertainties into account, it is far from clear whether the apparent excess crustal thickness beneath the Gulf is real; it is certainly close to the noise level in the data, and may thus instead be an artefact of the data processing rather than a real effect. It is unfortunate that Tiberi et al. (2001) did not publish their raw data, which would make it easier to test such possibilities, but until this issue is resolved (which is beyond the scope of this study) further gravity modelling of this region seems pointless. The nature of the technique used, combined with the assumption of a flat Moho, means that the numerical modelling solutions presented in this study are consistent with a difference in Bouguer gravity anomaly between the Gulf and its surroundings to the south of zero. The likely difference in behaviour, if the flat Moho constraint imposed on the modelling program for the solutions in Figs. 5, 6 and 7 were removed, can be inferred by analogy with solutions by Morley and Westaway (2006) and Westaway et al. (2006b) for vertical crustal motions, which persist over longer spatial scales, in other regions. The expected behaviour would be for the Moho to become deflected upward

81

beneath the Gulf depocentre and downward beneath the PGCB sediment source. For a given value of ηe, the flux of lower crust from beneath the depocentre to beneath the sediment source would be roughly the same as before, but the adjustment process occurring in the mantle would mean that this movement of material would result in less subsidence in the Gulf and less uplift in the PGCB. To match the observational evidence one would thus reduce the assumed value of ηe slightly, consistent with the assumption of a slightly higher Moho temperature. Such adjustments would not have any significant effect on the main conclusions, so no solutions of this type are presented in this study. As Morley and Westaway (2006) have noted, the question whether or not depocentres are underlain by Moho keels (i.e., underlain by the Moho at greater depths than beneath the surroundings to the basin) provides a significant observational test for lowercrustal flow driven by sediment loading. The presence of an “anti-keel”, where the Moho is shallower beneath the basin than beneath its surroundings, is consistent with lower-crustal flow in this sense, provided of course that the scale of the basin is large compared with the flexural wavelength of the underlying mantle lithosphere (otherwise the Moho will remain flat). Knowledge of the detailed Moho topography beneath the Gulf of Corinth is thus of some significance for testing the general class of model (Fig. 4) that is developed in this study. If the Moho beneath the Gulf is flat or an “antikeel” is present, then the proposed general form of solution is tenable, albeit in the latter case with minor modifications. Conversely, if a Moho keel could be shown to be present, then the concept would be invalidated. For reasons already discussed, I do not consider the evidence for the Moho “keel”, reported by Tiberi et al. (2001), to be compelling; this issue will only be resolved, following the development of much better data for crustal thickness variations, in the future. Another observational test may eventually be possible for the Gulf of Corinth, based on comparison of the ∼ 4 × 1019 Pa s estimate for ηe, from the present modelling, with estimates for the viscosity of the lower crust derived from geodetic study of post-seismic deformation following future large earthquakes. Such a test has already been carried out for western Turkey by Westaway et al. (2006b), whose numerical modelling of the Quaternary uplift history at Kula indicates ηe ∼ 2 × 1018 Pa s. Using the Westaway (1998) calibration, this implies a Moho viscosity (i.e., a viscosity in the deepest crust, just above the Moho) of ∼3 × 1016 Pa s and a Moho temperature of ∼660 °C in the Kula area. For comparison, geodetic studies of post-seismic deformation

82

R. Westaway / Tectonophysics 440 (2007) 67–84

following the August 1999 İzmit earthquake in NW Turkey (e.g., Hearn et al., 2002) indicate a Moho viscosity of ∼1017 Pa s, suggesting a slightly lower Moho temperature, consistent with the lower surface heat flow evident in NW Turkey relative to central western Turkey (cf. Tezcan, 1979, 1995; Ilkışık, 1995; Pfister and Rybach, 1996; Pfister et al., 1998). The fact that the lower-crustal flow occurring on time-scales that differ by so many orders-of-magnitude indicates consistent viscosities supports the view that the underlying rheology is linear, as previously assumed in my own investigations (following Westaway, 1998), rather than a power-law as is often assumed by others. The success of such comparisons is accompanied by a growing body of evidence that lower-crustal flow, induced by surface processes and thus driven by longtimescale climate change, can account for much of the growth of surface relief that has occurred worldwide in the Late Cenozoic (e.g., Bridgland and Westaway, in press). It is nonetheless evident that some people do not accept that lower-crustal flow is significant, within the continental lithosphere in general, except possibly in extreme localities such as Tibet where the Moho temperature is abnormally high (e.g., McKenzie et al., 2000; Maggi et al., 2001; McKenzie and Jackson, 2002). Jackson (2002) indeed claimed that the lower continental crust is typically so ‘stiff’ that it is ‘stronger’ than the mantle lithosphere, a view that is contrary to a wealth of evidence and not itself supported by any proper evidence. This sequence of publications has indeed been severely criticised by Morley and Westaway (2006), not least for the underlying element of circular reasoning; this criticism is not repeated here. A major issue leading to this difference of views concerns choice of the rheological law governing flow within the lower crust. The use of power-law rheologies for deformation at geological strain rates, based on simple extrapolation of high-strain-rate laboratory measurements, is strongly engrained within Earth Science. Such an extrapolation will result in an infinite viscosity in the limit of zero strain rate, and thus leads to the prediction of very high viscosities at geological strain rates. Westaway (1998) proposed instead a more sophisticated form of extrapolation from laboratory data, based on a method widely used in Materials Science to blend power-law behaviour at laboratory strain rates with linear behaviour at geological strain rates. This results in the prediction of much lower viscosities at the strain rates and temperatures appropriate for the lower crust than the alternative approach. As already noted, the growing ability of modelling studies that incorporate this rheological assumption to explain observational

evidence increasingly supports the conclusion that this rheological approach is correct. If, on the contrary, one were to deny that lower-crustal flow is feasible under the conditions that exist beneath the Gulf of Corinth, then clearly the present study would be considered invalid, as would its predecessor by Westaway (2002a). This would mean that the most recent valid reference concerned with explaining the magnitude of post-Early Pleistocene uplift south of the Gulf of Corinth would be that by Armijo et al. (1996). The flexural modelling by Armijo et al. (1996) required the ratio of footwall uplift south of the Gulf to hangingwall subsidence must be ∼ 1:3, requiring the Gulf to be underlain by as much as ∼ 5 km of sediment. However, the observed thickness of sediment beneath the Gulf was already known to be much less than this, thereby motivating the modelling by Westaway (2002a) to test if lower-crustal flow can explain the observed morphology instead. In the meantime, Clement et al. (2004) have confirmed that the maximum sediment thickness beneath the Gulf is no more than ∼ 1.7 km. Since the modelling by Armijo et al. (1996) can thereby be seen to be wrong, if one discounts solutions incorporating lower-crustal flow, the only options left to account for the observed morphology remain the various older ad hoc ideas mentioned earlier; otherwise the extent of post-Early Pleistocene uplift of the northern Peloponnese must still be regarded as a complete mystery. 5. Conclusions Numerical modelling, incorporating coupling between surface processes and induced flow in the lower continental crust, has been used to address the Quaternary evolution of the Gulf of Corinth region in central Greece. The post-Early Pleistocene marine depocentre beneath this Gulf overlies the northern margin of an older (Early Pleistocene and earlier) lacustrine basin, the Proto Gulf of Corinth Basin or PGCB. In the late Early Pleistocene, relief in this region was minimal but, subsequently, dramatic relief has developed, involving the creation of ∼ 900 m of bathymetry within the Gulf and the uplift by many hundreds of metres of the part of the PGCB, south of the modern Gulf, which forms the Gulf's main sediment supply. It is assumed that, as a result of climate change around 0.9 Ma, erosion of this sediment source region and re-deposition of this material within the Gulf began, both processes occurring at spatial average rates of ∼ 0.2 mm a− 1. Modelling of the resulting isostatic response indicates that the local effective viscosity of the lower crust is ∼ 4 × 1019 Pa s, indicating a Moho

R. Westaway / Tectonophysics 440 (2007) 67–84

temperature of ∼560 °C. It predicts that the ∼10 mm a− 1 of extension across this ∼ 70 km wide model region, at an extensional strain rate of ∼ 0.15 Ma − 1 , is partitioned with ∼3 mm a− 1 across the sediment source, ∼ 2 mm a− 1 across the depocentre, and ∼5 mm a− 1 across the ‘hinge zone’ in between, the latter value being an estimate of the extension rate on normal faults forming the major topographic escarpment at the southern margin of the Gulf. This modelling confirms the previous suggestion, by Westaway (2002a), that coupling between this depocentre and sediment source by lower-crustal flow can explain the dramatic development in local relief since the late Early Pleistocene. The effective viscosity of the lower crust in this region is not particularly low; the strong coupling interpreted between the sediment source and depocentre results instead from their close proximity. In detail, the effective viscosity of the lower crust is expected to decrease northward across this model region, due to the northward increase in exposure of the base of the continental lithosphere to the asthenosphere; in the south the two are separated by the subducting Hellenic slab. The isostatic consequences of such a lateral variation in viscosity provide a natural explanation for why, since ∼0.9 Ma, the modern Gulf has developed asymmetrically over the northern part of the PGCB, leaving the rest of the PGCB to act as its sediment source. Acknowledgments This work contributes to International Geoscience Programme 518 ‘Fluvial sequences as evidence for landscape and climatic evolution in the Late Cenozoic’. I thank Ioannis Koukouvelas and two anonymous reviewers for their thoughtful and constructive comments. References Allen, M., Jackson, J., Walker, R., 2004. Reply to comment by Rob Westaway on “Late Cenozoic reorganisation of the Arabia–Eurasia collision and the comparison of short-term and long-term deformation rates”. Tectonics 23. doi:10.1029/2004TC001695. Armijo, R., Meyer, B., King, G.C.P., Rigo, A., Papanastassiou, D., 1996. Quaternary evolution of the Corinth Rift and its implications for the Late Cenozoic evolution of the Aegean. Geophysical Journal International 126, 11–53. Bridgland, D., Westaway, R., in press. Climatically controlled river terrace staircases: a worldwide Quaternary phenomenon. Geomorphology. Bridgland, D.R., Philip, G., Westaway, R., White, M., 2003. A long Quaternary terrace sequence in the Orontes River valley, Syria: a record of uplift and of human occupation. Current Science (New Delhi) 84, 1080–1089.

83

Bridgland, D.R., Demir, T., Seyrek, A., Pringle, M., Westaway, R., Beck, A.R., Rowbotham, G., Yurtmen, S., in press. Dating Quaternary volcanism and incision by the River Tigris at Diyarbakır, SE Turkey. Journal of Quaternary Science. Clement, C., Sachpazi, M., Charvis, P., Graindorge, D., Laigle, M., Hirn, A., Zafiropoulos, G., 2004. Reflection–refraction seismics in the Gulf of Corinth; hints at deep structure and control of the deep marine basin. Tectonophysics 391, 97–108. Collier, R.E.L., 1990. Eustatic and tectonic controls upon Quaternary coastal sedimentation in the Corinth basin, Greece. Journal of the Geological Society (London) 147, 301–314. Demir, T., Westaway, R., Bridgland, D., Pringle, M., Yurtmen, S., Beck, A., Rowbotham, G., in press. Ar-Ar dating of Late Cenozoic basaltic volcanism in northern Syria: implications for the history of incision by the River Euphrates and uplift of the northern Arabian Platform. Tectonics. Demir, T., Yeşilnacar, İ., Westaway, R., 2004. River terrace sequences in Turkey: sources of evidence for lateral variations in regional uplift. Proceedings of the Geologists' Association 115, 289–311. Doutsos, T., Piper, D.J.W., 1990. Listric faulting, sedimentation, and morphological evolution of the Quaternary eastern Corinth rift, Greece: first stages of continental rifting. Geological Society of America Bulletin 102, 812–829. Fytikas, M.D., Kolios, N.P., 1979. Preliminary heat flow map of Greece. In: Cermak, V., Rybach, L. (Eds.), Terrestrial Heat Flow in Europe. Springer-Verlag, Berlin, pp. 197–205. Hearn, E.H., Bürgmann, R., Reilinger, R.E., 2002. Dynamics of İzmit earthquake postseismic deformation and loading of the Düzce earthquake hypocentre. Bulletin of the Seismological Society of America 92, 172–193. Jackson, J., 2002. Strength of the continental lithosphere; time to abandon the jelly sandwich? GSA Today 12 (9), 4–10. Ilkışık, O.M., 1995. Regional heat flow in western Anatolia using silica temperature estimates from thermal springs. Tectonophysics 244, 175–184. Keraudren, B., Sorel, D., 1987. The terraces of Corinth (Greece) — a detailed record of eustatic sea-level variations during the last 500,000 years. Marine Geology 77, 99–107. Lachenbruch, A.H., 1970. Crustal temperature and heat production: implications of the linear heat flow relation. Journal of Geophysical Research 75, 3291–3300. Laigle, M., Sachpazi, M., Hirn, A., 2004. Variation of seismic coupling with slab detachment and upper plate structure along the western Hellenic subduction zone. Tectonophysics 391, 85–95. Le Pichon, X., Angelier, J., 1981. The Aegean Sea. Philosophical Transactions of the Royal Society, London. Series A 300, 357–372. Le Pourhiet, L., Burov, E., Moretti, I., 2003. Initial crustal thickness geometry controls on the extension in a back arc domain: case of the Gulf of Corinth. Tectonics 22 (4), 1032. doi:10.1029/2002TC001433 (16 pp. (published online)). Maddy, D., Demir, T., Bridgland, D., Veldkamp, A., Stemerdink, C., van der Schriek, T., Westaway, R., 2005. An obliquity-controlled Early Pleistocene river terrace record from western Turkey? Quaternary Research 63, 339–346. Maggi, A., Jackson, J.A., McKenzie, D., Priestley, K., 2001. Earthquake focal depths, effective elastic thickness, and the strength of the continental lithosphere. Geology 28, 495–498. Mariolakos, I., Stiros, S.C., 1987. Quaternary deformation of the Isthmus and Gulf of Corinthos (Greece). Geology 15, 225–228. McClusky, S., et al., 2000. Global positioning system constraints on plate kinematics and dynamics in the eastern Mediterranean and Caucasus. Journal of Geophysical Research 105, 5695–5719.

84

R. Westaway / Tectonophysics 440 (2007) 67–84

McKenzie, D.P., 1978. Some remarks on the development of sedimentary basins. Earth and Planetary Science Letters 40, 25–32. McKenzie, D.P., Jackson, J., 2002. Conditions for flow in the continental crust. Tectonics 21 (6). doi:10.1029/2002TC001394 (200212) 7 pp. (published online). McKenzie, D.P., Nimmo, F., Jackson, J.A., Gans, P.B., Miller, E.L., 2000. Characteristics and consequences of flow in the lower crust. Journal of Geophysical Research 105, 11029–11046. McNeill, L.C., Collier, R.E.L., 2004. Uplift and slip rates of the eastern Eliki fault segment, Gulf of Corinth, Greece, inferred from Holocene and Pleistocene terraces. Journal of the Geological Society (London) 161, 81–92. Morley, C.K., Westaway, R., 2006. Subsidence in the super-deep Pattani and Malay basins of Southeast Asia: a coupled model incorporating lower-crustal flow in response to post-rift sedimentloading. Basin Research 18, 51–84. Mudelsee, M., Schulz, M., 1997. The Mid-Pleistocene climate transition: onset of 100 ka cycle lags ice volume build-up by 280 ka. Earth and Planetary Science Letters 151, 117–123. Ori, G.G., 1989. Geologic history of the extensional basin of the Gulf of Corinth (? Miocene–Pleistocene), Greece. Geology 17, 918–921. Pfister, M., Rybach, L., 1996. High-resolution temperature logging in shallow drillholes for the determination of terrestrial heat flow; field examples and analysis. Tectonophysics 257, 93–99. Pfister, M., Rybach, L., Simsek, S., 1998. Geothermal reconnaissance of the Marmara Sea region (NW Turkey); surface heat flow density in an area of active continental extension. Tectonophysics 291, 77–89. Rettenmaier, D., Giurgea, V., Förster, A., Hötzl, H., 2005. The AIG10 drilling project (Northern Peloponnesus, Greece): from geological investigations, geophysical logging, and hydrogeological testing to modelling, 3 pp. published online; available from: www.gfz-potsdam. de/pb5/pb52/projects/Corinth/Rettenmaier_ICDP_abstract2005.pdf. Stewart, I., Vita-Finzi, C., 1996. Coastal uplift on active normal faults: the Eliki Fault, Greece. Geophysical Research Letters 23, 1853–1856. Tezcan, A.K., 1979. Geothermal studies, their present status and contribution to heat flow contouring in Turkey. In: Cermak, V., Rybach, L. (Eds.), Terrestrial Heat Flow in Europe. SpringerVerlag, Berlin, pp. 283–292. Tezcan, A.K., 1995. Geothermal explorations and heat flow in Turkey. In: Gupta, M.L., Yamano, M. (Eds.), Terrestrial Heat Flow and Geothermal Energy in Asia. Science Publishers, Lebanon, New Hampshire, pp. 23–42. Tiberi, C., Lyon-Caen, H., Hatzfeld, D., Achauer, U., Karagianni, E., Kiratzi, A., Louvari, E., Panagiotopoulos, D., Kassaras, I., Kaviris, G., Makropoulos, K., Papadimitriou, P., 2000. Crustal and upper mantle structure beneath the Corinth Rift (Greece) from a teleseismic tomography study. Journal of Geophysical Research 105 (28), 28,159–28,171.

Tiberi, C., Diament, M., Lyon-Caen, H., King, T., 2001. Moho topography beneath the Corinth Rift area (Greece) from inversion of gravity data. Geophysical Journal International 145, 797–808. Westaway, R., 1996. Quaternary elevation change in the Gulf of Corinth of Central Greece. Philosophical Transactions of the Royal Society, London. Series A 354, 1125–1164. Westaway, R., 1998. Dependence of active normal fault dips on lowercrustal flow regimes. Journal of the Geological Society (London) 155, 223–253. Westaway, R., 2002a. The Quaternary evolution of the Gulf of Corinth, central Greece: coupling between surface processes and flow in the lower continental crust. Tectonophysics 348, 269–318. Westaway, R., 2002b. Geomorphological consequences of weak lower continental crust, and its significance for studies of uplift, landscape evolution, and the interpretation of river terrace sequences. Netherlands Journal of Geosciences 81, 283–304. Westaway, R., 2004. Comment on “Late Cenozoic reorganization of the Arabia–Eurasia collision and the comparison of short-term and long-term deformation rates” by M. Allen, J. Jackson, R. Walker. Tectonics 23, TC5006. doi:10.1029/2004TC001674 4 pp. (published online). Westaway, R., 2006a. Cenozoic cooling histories in the Menderes Massif, western Turkey, may be caused by erosion and flat subduction, not low-angle normal faulting. Tectonophysics 412, 1–25. Westaway, R., 2006b. Late Cenozoic extension in southwest Bulgaria: a synthesis. In: Robertson, A.H.F., Mountrakis, D. (Eds.), Tectonic Development of the Eastern Mediterranean Region. Geological Society, 260. Special Publications, London, pp. 557–590. Westaway, R., 2006c. Investigation of coupling between surface processes and induced flow in the lower continental crust as a cause of intraplate seismicity. Earth Surface Processes and Landforms 31, 1480–1509. Westaway, R., Bridgland, D., in press. Late Cenozoic uplift of southern Italy deduced from fluvial and marine sediments: coupling between surface processes and lower-crustal flow. Quaternary Intenational. Westaway, R., Pringle, M., Yurtmen, S., Demir, T., Bridgland, D., Rowbotham, G., Maddy, D., 2004. Pliocene and Quaternary regional uplift in western Turkey: the Gediz river terrace staircase and the volcanism at Kula. Tectonophysics 391, 121–169. Westaway, R., Demir, T., Seyrek, A., Beck, A., 2006a. Kinematics of active left-lateral faulting in southeast Turkey from offset Pleistocene river gorges: improved constraint on the rate and history of relative motion between the Turkish and Arabian plates. Journal of the Geological Society (London) 163, 149–164. Westaway, R., Guillou, H., Yurtmen, S., Beck, A., Bridgland, D., Demir, T., Scaillet, S., Rowbotham, G., 2006b. Late Cenozoic uplift of western Turkey: improved dating of the Kula Quaternary volcanic field and numerical modelling of the Gediz river terrace staircase. Global and Planetary Change 51, 131–171.