The structure of the lithosphere and upper mantle beneath the Eastern Mediterranean and Middle East

Surface velocity measurements show that the Middle East is one of the most actively deforming regions of the continents. The structure of the underlying lithosphere and convecting upper mantle can be explored by combining three types of measurement. The gravity field from satellite and surface measurements is supported by the elastic properties of the lithosphere and by the underlying mantle convection. Three dimensional shear wave velocities can be determined by tomographic inversion of surface wave velocities. The shear wave velocities of the mantle are principally controlled by temperature, rather than by composition. The mantle composition can be obtained from that of young magmas. Application of these three types of observation to the Eastern Mediterranean and Middle East shows that the lithosphere thickness in most parts is no more than 50-70 km, and that the elastic thickness is less than 5 km. Because the lithosphere is so thin and weak the pattern of the underlying convection is clearly visible in the topography and gravity, as well as controlling the volcanism. The convection pattern takes the form of spokes: lines of hot upwelling mantle, joining hubs where the upwelling is three dimensional. It is the same as that seen in high Rayleigh number laboratory and numerical experiments. The lithospheric thicknesses beneath the seafloor to the SW of the Hellenic Arc and beneath the NE part of the Arabian Shield are more than 150 km and the elastic thicknesses are 30–40 km.


Introduction
Since its discovery in 1967 it has been clear that plate tectonics only provides a good description of oceanic plate boundaries. On continents earthquakes are not confined to narrow belts (Isacks et al. 1968), and this distribution is not the result of location errors. This difference between oceans and continents is why plate tectonics was not discovered by geologists working on land. One of the first applications of the new methods to continental deformation was to the western part of the Alpine-Himalayan Belt by McKenzie (1972), who showed that most of the convergence between Africa and Europe in the Eastern Mediterranean was taken up by the rapid westward motion of Anatolia and by subduction of the Eastern Mediterranean beneath the Hellenic Arc, and not by shortening on east-west striking thrusts. This behaviour accounted for the strike slip motion on the North and East Anatolian Faults that had been mapped by Ketin (1948) and Şaroğlu (see Emre et al. 2018). As Ketin showed, the lateral motion on east-west striking North Anatolian Fault is at right angles to the slip vector between the major plates of Africa and Europe. Its motion was therefore not the expression of the relative motion between the large plates, and was not obviously compatible with the relative velocities obtained from the oceanic magnetic anomalies. Instead it results from the rapid westward motion of the Anatolian Plate, a small plate that is being driven westwards by boundary forces and gravitation. The eastern part of Tibet was later found to be behaving in the same way (Molnar and Tapponnier 1975), though it does not move eastwards as a rigid plate like Anatolia. Such motions must have been widespread in the past, and must have removed small pieces 1 3 of oceanic lithosphere that became trapped between colliding continents where their margins do not match each other.
Though the westward motion of Anatolia was first discovered from the surface displacements on active faults and the slip vectors of earthquakes, it has been spectacularly confined by GPS and InSAR measurements (e.g. Weiss et al. 2020). These have shown that the strain resulting from the motion of Anatolia is confined to a zone centred on the trace of the fault whose width is about 60 km, and that the interior of the Anatolian Plate is scarcely being deformed or elastically strained. In Greece the westward motion of Anatolia is partly taken up by bookshelf faulting, on many normal faults with a left lateral strike motion which bound rotating blocks (McKenzie 1978). D' Agostino et al. (2020) have used GPS measurements to show that the present rate of clockwise rotation of fault-bounded blocks in Greece is similar to that obtained from paleomagnetism.

Background
The surface horizontal velocities on Earth are now relatively well constrained, both on land and at sea, using magnetic anomalies and the strikes of fracture zones. But the resulting model is purely kinematic, and only provides a boundary condition, of the horizontal velocity at the surface, on the circulation of the mantle below, whose dynamics is controlled thermal buoyancy and by its mechanical properties. Beneath the oceans deformation is largely confined to narrow plate boundaries, and thermal models of the rigid cooling lithosphere account for much of the observed variation in water depth. Furthermore the elastic thickness, T e , of much of the oceanic lithosphere is 20-30 km, which is sufficiently large to suppress the surface expression of any mantle convection whose wavelength is less than about 500 km. These properties of the oceanic lithosphere explain why plate tectonics can be so successfully applied to the oceans. However, continents behave rather differently. Large regions have lithospheric thicknesses of less than 100 km and elastic thicknesses of no more than 5 km. Within such regions the surface elevation, the gravity field and the volcanism reflect any convective motions in underlying mantle. In contrast, many cratons have lithospheric thicknesses of 150-200 km and elastic thicknesses of 20-30 km, causing them to resist any tectonic deformation. Because of their lithospheric thickness, there is little surface expression of any convective circulation that is occurring in the underlying mantle.
The creep properties, or rheology, of solid materials are well understood by those who study the properties of materials (see Stocker and Ashby 1973). They depend on the homologous temperature ∕ s and the homologous stress ∕ , where is the temperature in Kelvin, s ≃ 1200 • C the solidus temperature, the applied shear stress and the shear modulus whose value in the upper mantle is about 7 × 10 10 Pa. In most places the horizontal motions at the Earth's surface extend into the mantle and the plate motions are taken up at shallow depths by brittle deformation. However, if the crust is thick and its lower part is hot, as it is beneath Tibet, the mantle beneath the Moho may be decoupled from the surface motions. At sufficiently high temperatures elastic stresses are relaxed by creep and the material behaves as a Maxwell solid rather than as an elastic solid with a yield stress. The time scale for stress relaxation m is then governed by m = ∕ , where is the viscosity. If the stresses involved in the generation of earthquakes are ∼ 20 MPa, and are produced on time scales of ∼ 10 3 a, they will be relaxed by creep if < m ≃ 2 × 10 21 Pa s. Stocker & Ashby's deformation map shows that this value of corresponds to a value of ∕ s ≃ 0.6 when ∼ 20 MPa. The corresponding values for = 5 MPa are ∕ s ≃ 0.65 and T = 680 • C. Therefore seismicity will be restricted to regions whose temperature is less than about 600 • C-700 • C. This simple argument gives values for the limiting temperature that is in the same range as that obtained from the depths of oceanic earthquakes (see McKenzie et al. 2019a). The stresses that support topography and flexure must be maintained for times of at least 1 Ma, or much longer than those involved in earthquakes, and involve stresses of a similar magnitude. The same approach then gives an estimate of the limiting viscosity of ∼ 2 × 10 24 Pa s and an upper limit to the temperature at which such stresses can be maintained of about 500 • C. Though these temperature estimates are not likely to be accurate, the difference in the time scales requires that the thickness, T e , of the layer that provides elastic support for flexure and topography must be less than that of the seimically active layer T s .
At temperatures greater than about 600 • C elastic stress are relaxed by creep, but at low stresses the viscosity is independent of the stress. For a value of ∕ s ≃ 0.60 ∼ 10 25 Pa s and m ∼ 4.5 Ma whereas ∕ s ≃ 0.65 gives ∼ 10 24 Pa s and m ∼ 0.5 Ma (Stocker and Ashby 1973). These estimates of the viscosity are sufficiently large that the material will move with the surface velocity. In most places this region includes the crust and upper mantle to a depth of at least 50-70 km, and is referred to as the mechanical boundary layer, or MBL. Within the MBL heat is transported by conduction. The time taken for the MBL to reach steady state depends on its thermal time constant, which in turn depends on the square of its thickness. Beneath the oceans the thickness of the MBL is about 100 km and its time constant is about 60 Ma. Where the MBL is as thick as 180 km, as it is beneath some cratons, the time constant reaches 200 Ma. In the region below the MBL, where the viscosity is sufficiently low to allow the mantle to move with respect to the surface, heat is transported by advection as well as by conduction. This region forms the thermal boundary layer, or TBL, to the convective circulation below. At the upper boundary of the TBL all heat is transported by conduction, whereas at its lower boundary heat transport is dominated by advection. Both the temperature and the temperature gradient must be continuous everywhere. For this reason none of the boundaries of the MBL and TBL are sharp. Though heat transport is dominated by advection in the interior of the convective circulation, conductive transport must always occur where there are temperature gradients. In the interior the horizontally averaged entropy is almost constant. It is therefore convenient to use the potential, rather than the actual, temperature in discussions of mantle temperature. The potential temperature is defined as the temperature the material would have if it were brought to the surface at constant entropy without melting.
From the geological point of view an important feature of mantle convection is the pattern of the convective circulation when viewed from above, which is known as its planform. Wilson (1963) and Morgan (1971) suggested that the planform of mantle convection consists of axisymmetric plumes and their proposal has been widely accepted. However, it was based on an intuitive approach, guided by geological and geophysical observations, rather than on laboratory or numerical studies of high Rayleigh number, high Prandtl number convection. Such laboratory experiments Whitehead 1971, 1974;Richter and Parsons 1975) show that the planform is more complicated than Wilson and Morgan supposed, and consists of localised upwellings, or hubs, joined together by two-dimensional rising sheets, or spokes, to form a planform that is known as a spoke pattern. It is the presence of hubs that led to Wilson and Morgan's proposals. Recent numerical experiments by Lees et al. (2020) in large aspect ratio boxes show that the spokes, as well as the hubs, can deform the surface, generate gravity anomalies and magma, provided the overlying lithosphere is sufficiently thin.
Clearly this framework, which depends on ideas from fluid mechanics and material science, should form the basis for the interpretation of geological and geophysical observations. An important difficulty in doing so concerns the concept of 'lithosphere' and 'lithosphere-asthenosphere boundary', or LAB. The definition of lithospheric thickness used in this paper is the depth where the conductive geotherm reaches a potential temperature of 1300 • C, which is at present the best estimate of the globally averaged potential temperature of the convecting upper mantle. This depth does not correspond to any physical interface within the mantle, and is in the middle of the thermal boundary layer. Therefore the thickness of the lithosphere defined in this way can obviously not be determined by any direct geophysical measurement, such as the depth at which the conversion of P to S waves occurs. A variety of other definitions of lithospheric thickness have been proposed in the half century since plate tectonics was discovered. Like that used here, none of them correspond to Barrell's definition (see Watts 2001), which would place the base of the lithosphere within the crust in regions like Tibet. The relationship between these other definitions and the thermal, mechanical and fluid dynamical model outlined above is not yet understood and is the subject of much present research.
Three types of measurement can be used to study the structure of the lithosphere and of the mantle circulation below: of gravity, of seismic velocities and of the compositions of magmas generated in the underlying mantle. Gravity measurements have long been used for this purpose (e.g. Heiskanen and Vening Meinesz 1958). Covering large areas using surface measurements is both slow and expensive. However, unless the continent is covered with water, surface measurement is still the only available method of determining anomalies with wavelengths of less than about 200 km that are required to estimate values of T e that are less than 5 km. But wherever the surface of the continent or oceanic crust is covered by seawater, measurement of the shape of the sea surface by altimetric satellites, which in most places is a good approximation to the geoid, now provides the most accurate and fastest method of obtaining the gravity field (e.g. Sandwell et al. 2013).
In subaerial regions gravity anomalies with wavelength of 200 km or more have now been accurately determined from satellite measurements, especially those from GRACE and GOCE (Bruinsma et al. 2014). The technology involved in making such measurements is remarkable. The GOCE satellite contained small masses separated by about 0.5 m, and measured the force that was required to maintain them at the centre of their cavities. For a gravity anomaly with a wavelength of 200 km and an amplitude of 100 mGal at the Earth's surface, the difference in the gravitational acceleration between such masses in a satellite at a height of 235 km is about 10 −12 g, where g is the acceleration due to gravity at the Earth's surface. Recent measurements of the relative velocity between two satellites using an optical laser (Ghobadi-Far et al. 2020) suggest that this system may be able to do even better than GOCE.
Two methods of estimating T e , the elastic thickness of the lithosphere, from gravity anomalies are used below. The first fits stacked profiles of free air gravity anomalies to the calculated profile of a plate flexed by the application of a force and a bending moment (see Heiskanen and Vening Meinesz 1958 ch. 10), using the approach of McKenzie et al. (2019b). This method works well where the lithosphere is loaded by a long thrust fault, but is difficult to use in regions of rough topography. In such regions T e is instead estimated using spectral methods, by finding the transfer function, Z(k), where k is the wavenumber, often referred to as the admittance, between the two dimensional Fourier transform of the topography as input and that of free air gravity as output (see Watts 2001 ch. 5). At long wavelengths, where flexural effects become unimportant, the admittance is often controlled by convection, and is about 50 mGals/km over land and 30 mGals/km over the sea. At wavelengths between 500 and 100 km the topography changes from being compensated to uncompensated. As it does so the behaviour of Z(k) can be used to estimate T e (McKenzie and Fairhead 1997).
Determination of seismic velocities depends on the measurement of ground motion by seismometers, which is then used to model the three dimensional velocity structure. The largest amplitude signal is that of the fundamental mode Rayleigh waves, whose amplitude at a distance of 3000 km from a magnitude 6 earthquake is between 10 and 100 μ m at a period of 50 s. Ground motion can now be accurately measured using seismometers that determine the force required to keep a proof mass stationary with respect to the ground, and GPS provides accurate timing. Conversion of such measurements into velocity models is an active area of research. Two approaches are widely applied, both of which can be used to generate regional as well as global models of V p and V s . One is based on the classical approach to seismology, and uses the measurement of the velocities of P, S, Rayleigh and Love waves between sources and a global network of seismometers. It is especially important to include measurements of the higher modes of Rayleigh and Love waves, because these are important in constraining the depth of any velocity anomalies. The velocities are then inverted to produce an isotropic or an anisotropic velocity model (Priestley et al. 2019;Schaeffer and Lebedev 2013) Such studies typically use 10 6 -10 7 seismograms and produce global models with a horizontal resolution of 200-300 km and a vertical resolution in the upper mantle of about 30 km. The agreement between different models, produced by different approaches to inversion, is now quite good (see supplementary material). This approach has the advantage of being able to use all available seismograms, and of being additive: each seismogram need only be processed once, and does not need to be reprocessed for each iteration. Its disadvantage is that it can only use those parts of the seismograms that are dominated by waves which travel along the great circle from the source to the receiver. It cannot use those parts of seismograms where multipathing has occurred. Multipathing is frequency dependent, and often dominates Rayleigh and Love waves with periods shorter than 30 s. The second approach, known as full waveform inversion, generates synthetic seismograms by solving the elastodynamic equations directly, often carried out using spectral element codes, which are directly compared with observed seismograms. Because it makes no distinction between different types of seismic wave, and solves the full governing equations without making approximations, it is untroubled by such problems as multipathing and P to S conversion. In principle this approach, which uses every wiggle on the seismograms, should be more powerful. But in practice it has important limitations (see Lei et al. 2020).
To produce global models with useful resolution requires the use of the most powerful computers that exist, and even then can only use a relatively small number of seismograms. Lei et al. (2020) used 1480 earthquakes, or three orders of magnitude fewer than Priestley et al. (2019) and Schaeffer and Lebedev (2013) used. Though the rapid improvement in computer power shows no sign of coming to an end, it is likely to continue to be overwhelmed by the continuing rapid increase in the number of digital seismometers generating seismograms. Therefore the spatial resolution of full waveform inversion is likely to remain less than that of the first approach. Another problem is that it is not additive: each iteration requires the recalculation of every seismogram that is used. Another serious limitation concerns the criteria used to match the synthetic to the observed seismograms. The largest amplitude signal on many seismograms is that of the fundamental mode Rayleigh wave. The amplitudes of higher modes, whose behaviour provides the major constraints on the velocity variation with depth, are typically an order of magnitude smaller than that of the fundamental mode. The behaviour of the higher modes is especially important, because in a half space they are only present when the velocity varies with depth. The travel times of the higher modes, but not of the fundamental mode, are similar, but their velocities respond in different ways to the vertical velocity structure. Therefore, if the information in the higher modes is to constrain the velocity model, a weighting must be applied to the misfit between the observed and calculated seismograms. Clearly the choice of the weighting will have a major influence on the resulting model, but there is no general theory which can be used choose a weighting function with the desired properties. The discussion below is based on a shear wave velocity model produced by Priestley et al. (2019) using the first approach.
However they are produced, isotropic models from seismology provide three dimensional models of the density, V p and V s . The resulting variations in V s are larger and better constrained than are those of density or V p . Laboratory experiments show that variations in V s are more sensitive to temperature than to likely compositional variations, especially at temperatures above 800 • C (Schutt and Lesher 2006). Priestley et al.'s (2019) global V s model was therefore converted to one of temperature, which was then used to obtain a map of lithospheric thickness (Priestley and McKenzie 2006) for comparison with the free air gravity maps and melt composition.
The third type of measurement that provides important constraints on the properties of the upper mantle is the composition of mafic magmas. Modern inductively coupled plasma mass spectrometers (ICP-MS) and X-ray flourescence instruments (XRF) can measure the concentrations of 40 or more elements, many of which are present only in ppm, to a relative accuracy and precision of 5% or better. Concentration measurements of all the rare earth elements (REE) are now usually reported and are easier to interpret than are those of most elements, principally because all have a valency of three in most geological systems. Their geochemical behaviour is controlled by the uniform reduction in ionic radius that occurs with increasing atomic number. Such elements are not in fact rare: they are present in all mafic rocks in concentrations of at least a few ppm. Without any measurements of composition, the presence of magma alone constrains the temperature at some depth to be above the solidus. The REE compositions can often be used to constrain the depth from which the magmas originate, especially when combined with geotherms from seismic tomography. Furthermore, unlike the seismic velocities, the variation of the REE concentrations with atomic number constrains the composition and density of the source material of the magmas. McKenzie and O'Nions's (1991) method, modified to use Blundy and Wood's (1994) parameterisation of partition coefficients of the REE, is used below to determine the melting conditions and the density of their source regions.
The discussion below depends on only these three types of measurement to constrain the structure of the lithosphere and upper mantle. It does not make use of two other types of observations that are commonly used for this purpose. The first are estimates of heat flow, from temperature gradients in boreholes. Such measurements are difficult to interpret, partly because such temperatures are strongly affected by any fluid movement, of water, oil or gas (and such holes are Fig. 1 a Free air gravity anomalies from Sandwell et al. (2013) from the Eastern Mediterranean. The profiles inside the magenta box were used to construct the profile in Fig. 3a and to estimate the value of T e of 29 km. b locations of regions discussed in the text often drilled for hydrological purposes or petroleum exploration) and partly because continental heat flow is often dominated by upper crustal radioactivity (Jaupart and Mareschal 2005). The second are estimates of the electrical conductivity. Except near or above the solidus, the bulk electrical conductivity of most minerals present in the crust and upper mantle is very low. Increased conductivity can result from a variety of physical processes, especially ones that affect ionic movements on grain boundaries. Therefore the interpretation of conductivity anomalies is not straight forward.
The next section, Sect. 3, is concerned with constraints that gravity measurements and models of shear wave velocity impose on the structure of the lithosphere beneath the Eastern Mediterranean and Middle East. Section 4 discusses the structure of the upper mantle beneath the lithosphere using the same material, together with the geochemical measurements from mafic magmas. Section 5 is concerned with the structure beneath a limited region of the NW Zagros, which resembles that beneath northern India and Tibet.
3 Lithospheric structure Figure 1 shows the free air gravity anomalies over the Eastern Mediterranean and surrounding regions from Sandwell et al. (2013). Where the surface is covered by seawater the gravity field is obtained from altimetric satellite measurements and is well determined at wavelengths longer than 100 km. On land the gravity field at wavelengths between 100 and 200 km depends on surface measurements, which are sparse in some regions of Fig. 1. The free air gravity anomalies are negative over most of the Eastern Mediterranean Sea, partly because it has been overriden by thrusts along several of its boundaries. Those associated with the Hellenic Arc to the north are now active. To the south the deep trough along the Libyan continental margin, Fig. 2, suggests that the thrusting there may still be active, though the seismicity is less than that of the Hellenic Arc. Gravity profiles from a box SW of Greece, marked by the magenta boxes in Figs. 1 and 2, were used to produce the stacked profile in Fig. 3a (McKenzie et al. 2019b). The value of T e that best fits the profile is 29 km, calculated using a density contrast between the mantle and overlying sediments of 900 kg m −3 . Though the misfit in Fig. 3b shows that the upper bound on T e is not well determined, the value of T e must be greater than about 20 km. The model for the profile obtained from the free air gravity east of Honshu, Fig. 3c, uses a density contrast of 2300 kg m −3 between the mantle and the overlying material, because, unlike the Hellenic Arc, the trench east of Honshu is filled with water, not sediment. The flexure SW of the Hellenic Arc in Fig. 3a is gentler and the sea floor is not offset by large normal faults as it is east of Honshu. The gravity anomalies to the SE of the eastern part of the Hellenic Arc are similar to those used to produce Fig. 3a. However, the stacked profile from this part of the Arc is less straightforward to model, because it is complicated by the positive anomaly associated with the load from the Nile Delta (Fig. 1). Two earthquakes have occurred beneath the seafloor south of the Hellenic Arc with depths of 33 and 38 km (Jackson personal communication), the first of which is inside the box used to estimate T e . Close to the Arc the depths increase to 54 km (Shaw and Jackson 2010). As expected from the discussion above, the depths of the earthquakes are greater than the value of T e estimated from flexure.
The lithosphere beneath the part of the Mediterranean marked by the magenta box in Figs. 1 and 2 is generally believed to be oceanic lithosphere at least Cretaceous or Jurassic in age, and possibly even as old as the late Palaeozoic (Granot 2016). Therefore it is sufficiently old to have reached a thermal steady state. It is overlain by at least 10 km of sediment. Figure 4 shows how its thermal evolution is affected by this sediment blanket. The seismogenic thickness of oceanic lithosphere is approximately controlled by the depth of the 600 • C isotherm (see McKenzie et al. 2019a). Figure 4a, calculated using the approach of McKenzie and Jackson (2020), shows that the upper part of the oceanic crust and upper mantle beneath the sediments remains sufficiently cool to generate earthquakes and support the flexural stresses even after 120 Ma. Subduction of this lithosphere beneath the Hellenic Arc generates a sinking slab, with earthquakes reaching a depth of 150 km (see Howell et al. 2017). In regions where the temperature of subducting slabs can be estimated, their seismicity is restricted to regions cooler than about 600 • C. Figure 4b shows that the lithosphere in Fig. 4a will remain cooler than 600 • C to a depth of 400 km when it is subducted at 40 mm a −1 . This depth is greater than that of the deepest events beneath the Aegean, probably because the subduction has not yet been occurring for long enough to have reached steady state. Figure 5 shows that a flexural anomaly similar to that in the Eastern Mediterranean is present to the SW of the Zagros. The stacked profile in Fig. 3e is best fit with a value of T e of 37 km, or the same as that from McKenzie et al. (2019b) for northern India shown in Fig. 3g. Figure 5 shows three regions of rough topography whose values of T e , of between 3.4 and 4.2 km, were estimated from values of the admittance (Fig. 6). All three regions contain active volcanoes whose magmas are generated at shallow depths (see Sect. 4). Because the values of T e are so small, the topography and gravity in the Arabian and Turkish boxes with wavelengths greater than about 300 km are dominated by convective signals, with admittances of about 50 mGals km −1 (Fig. 6). Figure 7 shows gravity anomalies from GOCE (Bruinsma et al. 2014) whose wavelengths are larger than 300 km. As Fig. 6c illustrates, where the elastic thickness is less than 5 km such anomalies result from convection beneath the lithosphere. Active volcanoes and young lava flows in Fig. 7 are shown in black, and most coincide with positive gravity anomalies that mark convective upwellings. Wilson (1963) and Morgan (1971) suggested that the planform of mantle convection consists of axisymmetric plumes and their proposal has been widely accepted. Various modifications (Ebinger and Sleep 1998;Sleep 1997Sleep , 2008 have been proposed to Wilson and Morgan's ideas to explain how linear features like those in Fig. 7 can be produced by axisymmetric plumes. However, the laboratory and numerical experiments discussed above show that the expected planform for high Rayleigh number convection is a spoke pattern, and that the spokes, as well as the hubs, can generate gravity Fig. 4 a Shows the evolution of the temperature structure of oceanic lithosphere, which has reached a steady state temperature distribution at the start of sedimentation at 0 Ma, when 10 km of sediment is deposited on it at a constant rate over 120 Ma. b shows the steady state temperature structure in a sinking slab when the lithosphere in a is subducted at 40 mm a −1 120 Ma after sedimentation started (modified from McKenzie and Jackson 2020) anomalies and magma provided the overlying lithosphere is sufficiently thin. However, unlike the hubs, the spokes are hot rising sheets, and can therefore produce surface features like those in eastern Arabia, Syria and Anatolia in Fig. 7. Figure 8 shows that much of the young volcanism in the Middle East and NE Africa has occurred where the lithosphere is thin and the elastic thickness is small. The gravity anomalies show that this volcanism results from upwelling along a spoke that joins the hubs beneath the Afar to eastern Anatolia and NW Iran, and another spoke that trends EW beneath Turkey and NW Iran. An exception to this pattern is the region centred on the point labelled 'ham' in Fig. 8, which is underlain by lithosphere whose thickness reaches 190 km. This region is discussed in the next section. Figure 8 shows that three extensive regions of thick lithosphere, beneath the Eastern Mediterranean, the Zagros and the Caspian Sea, are all regions of negative gravity anomalies. Lithospheric thicknesses calculated from Schaeffer and Lebedev (2013) V s model, shown in the supplementary material, show the same features in the same locations, though the lithospheric thicknesses from their velocity model are somewhat thicker than those in Fig. 8. At least two of these locations are also regions where T e ∼ 30 km. Figure 6c shows that such large values of T e can provide elastic support for gravity anomalies with wavelengths of 400-500 km. It is therefore not clear whether the gravity anomalies in the Eastern Mediterranean are supported by flexure or, as Howell et al. (2017) proposed, by convection. In contrast, the positive gravity anomalies that extend across Turkey, through Syria, along the mountains north-east of the Red Sea to the Afar, and which are also associated with Darfur and Tibesti, where T e is everywhere less than about 5 km, are all likely to be convectively supported.

Convective structures below the lithosphere
There have been a number of local studies of seismic wave propagation in the region illustrated in Fig. 8. Blom et al. (2020) used full waveform inversion to obtain three dimensional models of V s , V p and for the central and eastern Mediterranean. In the regions where their resolution is good their results agree well with the global models and have better resolution than those used to produce Fig. 8. However, though they find some indications of its presence, their ray coverage is insufficient to map the large region of thick lithosphere that Priestley et al. (2019), Schaeffer and Lebedev (2013) and El-Sharkawy et al. (2020) find south of the Hellenic Arc. This region of thick lithosphere is the most enigmatic structure in the whole of the Mediterranean and Middle East. It is present in all tomographic models from surface waves. No similar feature, with lithospheric thickness of ∼ 170 km benath ∼ 4 km of water, has yet been discovered elsewhere on Earth. In the early low resolution tomographic models this anomaly was thought to result from the high velocities within the cold sinking slab beneath the Aegean, even though its location was south of the Hellenic Arc. However, as the resolution of the tomographic models has improved, it has become clear that this feature is indeed south of the Arc. Figure 8 shows that it has a similar structure to cratonic lithosphere, like that beneath SE Arabia, and does not resemble either that of the slab beneath NE Honshu, or that of the Mesozoic lithosphere of the NE Pacific. Figure 7 shows that the region has a large regional Fig. 6 Admittances and estimates of T e for the boxes in Figs. 5 and 10. That expected for convective support, of 50 mGals km −1 , is shown as a dotted line. The misfits as functions of T e are shown in the boxes on the right ◂ Fig. 7 Gravity from DIR-R5, smoothed to 300 km. The cyan points show the locations of the magma suites whose REE concentrations are modeled to obtain the melting conditions in Fig. 9, and the numbers in boxes show the best-fitting potential temperatures to the melt models in • C. gravity anomaly, of ∼ 120 mGals, superimposed on the flexural anomalies associated with the Hellenic Arc and the Libyan margin (Fig. 1). If the anomaly in Fig. 7 results from convection, as Howell et al. (2017) proposed, the corresponding surface depression should be about 4 km, in agreement with that observed in Fig. 2. Therefore the structure appears to be a piece of thick cratonic lithosphere that has been pulled down 4 km below sea level, either by local convective downwelling as Woodside (1976) suggested, or by the negative convective buoyancy transmitted through the lithosphere from the slab beneath the Aegean. Further east Skobeltsyn et al. (2014) and Kaviani et al. (2020) have used the dispersion of fundamental-mode Rayleigh waves, either from earthquakes or from noise, to produce three dimensional models of the V s structure beneath much of the Middle East. They found low velocities in regions where the lithosphere in Fig. 8 is thin and where there is extensive volcanism. The profiles labelled BB ′ and CC ′ in Figs. 9 and 10 of Skobeltsyn et al. (2014) probably cross part of the high velocity region corresponding to the thick lithosphere in the NW Zagros in Fig. 8.
Where the lithosphere and elastic thickness are large, convection is still likely to be present but with strongly reduced convective expression of the gravity and topography. McKenzie's (2010) expressions show that gravitational and topographic anomalies of 50 mGal and 1 km with a wavelength of 400 km, will be reduced to amplitudes of about 12 mGal and 200 m if T e = 30 km and the lithospheric thickness is 170 km. Anomalies resulting from flexure will then dominate those from convection, and the absence of linear positive gravity anomalies may simply result from the thickness of the lithosphere, rather than from the absence of spoke pattern convection. Figure 9 shows the REE concentration ratios from four suites of basalt from the locations in Fig. 7 marked by cyan circles, together with the melting models, and some of the geotherms calculated from surface wave tomography that were used to obtain the lithospheric thickness in Fig. 8. The models show that all the magmas can be generated by modest amounts of melting in the transition zone from garnet to spinel peridotite, in the regions where the geotherms exceed the solidus. The geotherms were calculated from the V s tomography, and the solidi from Katz et al.'s (2003) parameterisation with 0.02% water in the source. The same parameterisation can be used to estimate the potential temperature of the upwelling mantle, using the melt fractions obtained from the inversion and shown in the panels (b), (e), (h) and (k) of Fig. 9. The values are plotted in Fig. 7, and show that the potential temperature of the spoke is greater in central Arabia and Syria than it is in Yemen. All the temperatures from the melting models exceed the best estimate of the globally averaged potential temperature of the upper mantle, of 1300 • C, which is obtained from the condition that upwelling mantle beneath fast spreading ridges should generate 7 km of crust.
The composition of the magmas in the Middle East associated with convective upwellings is similar to that of ocean island basalts, and has a different composition from MORB. The upper mantle beneath oceanic lithosphere must therefore be inhomogeneous. The composition of the magmas in Fig. 8 is likely to result from the melting of inhomogeneities that are more fusible than is the bulk rock. The presence of such inhomogeneities can also account for the large variations in the volume of the eruptions shown in Fig. 7.
As well as providing estimates of the potential temperature, the REE composition of magmas can be used to estimate the modal composition of the sources of the magmas, and hence of their density. Figure 8 shows the difference between the density s of the sources of many of the magmas in the Middle East and that a of the MORB source at the same temperature, where = s − a and the units are kg m −3 . Except for the value at the point labelled 'ham', which is discussed in the next section, the calculated values show that the sources are slightly heavier than the MORB source, because, like the sources of ocean island basalts, they are slightly more fertile. If the sources were homogeneous they would need to be about 80 • C hotter than the MORB source to have the same density. In practice the sources are unlikely to be homogeneous and the magmas probably come from the more fusible parts of the mantle.

NW Zagros
Most of the seismic activity of the Zagros consists of moderately sized thrust earthquakes, with magnitudes of no more than about 6, which occur near the base of the sediments overlying the basement. The sediment thickness is about 12 km and the dip of most of the thrusts is about 45 • . These events are associated with the NE-SW shortening that produces the enormous folds. However the structure of part of the northern region illustrated in Fig. 10 closely resembles that of northern India and southern Tibet (see Craig et al. 2012), though it is on a smaller scale. In 2017 a magnitude 7.3 earthquake occurred in the region at a depth of 15 km on a large almost flat thrust (Feng et al. 2018;Yang et al. 2018, Fig. 10). The gravity field in Fig. 5 shows that the lithosphere to the SW of this earthquake has been flexed Fig. 9 Models of melt generation from modeling REE concentrations, and geotherms from surface wave velocity inversions. The black dots in c, f, i, l show the temperature estimates from surface wave tomography. The part of the geotherms marked by thick cyan lines shows the region where the melt is generated. a Uses measurements from Karacadaǧ from Ekici et al. (2014), d from South Syria from Krienitz et al. (2007), g from Harrat Khaybar from Camp et al. (1991), j from Yemen from Baker et al. (1997) ◂ downwards, and modeling of the stacked profiles in Fig. 3e gives a value of T e of more than 20 km, or similar to that of the flexure of northern India shown in Fig. 3g. Figure 8 shows that the region to the NE of the location of the earthquake is underlain by thick lithosphere, even though the value of T e from the gravity and topography is only 4 km (Fig. 6e, f). This small value suggests that, like Tibet, the hot lower crust decouples the upper crust from the thick underlying lithosphere. A further resemblance to the structure of northern India and Tibet is the composition of young volcanics from the region marked 'ham' in Fig. 10 (Allen et al. 2013). Magmas in both regions come from mantle material that was first depleted by the removal of 20-25% melt to form a harzburgite, which was then enriched by the addition of a few % of metasomatic melt, which was itself produced by the extraction of 1% or less of melt from the upper mantle. In both regions the source density is about 60 kg m −3 lighter than the average upper mantle at the same temperature (McKenzie and Priestley 2016). Though this density difference is not large, it is sufficient to cause the cold lithosphere to be less dense than the convecting upper mantle if its temperature is at least 800 • C, when it is convectively stable.
The composition of a magma from near Kerman, at 29.4 • N, 57.4 • E, is also similar to that in Fig. 10 and its source has a similar density deficit with respect to the average upper mantle (McKenzie and Priestley 2016) to that in Fig. 10 from near Hamadan. It is also underlain by thick lithosphere.

Conclusions
Geodetic, seismological and geochemical data obtained from new instruments, together with major advances in computer power, can be combined to model the structure of the lithosphere and the convecting mantle below. In the Eastern Mediterranean and Middle East most of the mountainous regions have small values of elastic thickness. Their topography with wavelengths of more than 300 km results from mantle convection in the form of lines of hot rising material, spokes, which link hubs (plumes). The upwelling associated with both spokes and hubs generate magmas at depths of 60-80 km. The melt fractions involved are less than 10% and the magma compositions are similar to those of ocean island basalts. The composition of the magmas is therefore controlled by that of the convecting upper mantle, and not by that of the overlying continental lithosphere. Fig. 10 The topography used to obtain estimates of T e is from Smith and Sandwell (1997), the gravity is from Eigen6c (Förste et al. 2011). The location of the magma whose composition was used to estimate the compositional density contrast is shown by a black dot, marked with in kg m −3 between its source and the average density of the upper mantle (see Fig. 8). The magma composition is from Allen et al. (2013). The focal mechanism is that of the Sarpol Zahab magnitude 7.3 earthquake on 2017.11.12 South of the Hellenic Arc and SW of the Zagros the elastic thickness of the lithosphere is greater than 20 km, and is probably as much as 30-40 km. This value is considerably greater than that of most continental regions which are not underlain by thick lithosphere. In these regions there is no obvious surface expression of any mantle convection, probably because the lithosphere is thick and strong rather than because such convection is absent.
A region in the NW Zagros southwest of the Caspian Sea resembles the structure of northern India and southern Tibet, though it is on a smaller scale. In this region the Zagros mountains are overriding and depressing the Arabian lithosphere on a large shallow angle thrust. Like Tibet, the overiding material has a small elastic thickness and is underlain by thick harzburgitic lithosphere. The presence of these features in both regions suggest that they result from extensive shortening of continental lithosphere, and may therefore have been commoner in the past than they are at present, since there are few continental regions now undergoing extensive shortening.