The Thickness of the Mantle Lithosphere and Collision-Related Volcanism in the Lesser Caucasus

The Lesser Caucasus mountains sit on a transition within the Arabia–Eurasia collision zone between very thin lithosphere ( < 100km) to the west, under Eastern Anatolia, and a very thick lithospheric root (up to 200km) in the east, under western Iran. A transect of volcanic highlands running from NW to SE in the Lesser Caucasus allows us to look at the effects of lithosphere thickness variations on the geochemistry of volcanic rocks in this continental collision zone. Volcanic rocks from across the region show a wide compositional range from basanites to rhyolites, and have arc-like geochemical characteristics, typiﬁed by ubiquitous negative Nb–Ta anomalies. Magmatic rocks from the SE, where the lithosphere is thought to be thicker, are more enriched in incompatible trace elements, especially the light rare earth elements, Sr and P. They also have more radiogenic 87 Sr/ 86 Sr, and less radiogenic 143 Nd/ 144 Nd. Across the region, there is no correlation between SiO 2 content and Sr–Nd isotope ratios, revealing a lack of crustal contamination. Instead, ‘spiky’ mid-ocean ridge basalt normalized trace element patterns are the result of derivation from a subduction-modiﬁed mantle source, which probably inherited its subduction component from subduction of the Tethys Ocean prior to the onset of continent–continent collision in the late Miocene. In more isotopically enriched mantle source, modelling of non-modal batch melting suggests lower degrees of melting and the involvement of garnet as a residual phase in the SE. Melt thermobarometry calculations based on bulk-rock major elements conﬁrm that melting in the SE must occur at greater depths in the mantle. Temperatures of melting below 1200 (cid:2) C, along with the subduction-modiﬁed source, suggest that melting occurred within the lithosphere. It is pro-posed that in the northern Lesser Caucasus this melting occurs close to the base of the very thin lithosphere (at a depth of (cid:3) 45km) as a result of small-scale delamination. A striking similarity between the conditions of melting in NW Iran and the southern Lesser Caucasus (two regions between which the difference in lithosphere thickness is (cid:3) 100km) suggests a common mechanism of melt generation in the mid-lithosphere ( (cid:3) 75km). The southern Lesser Caucasus magmas result from mixing between partial melts of deep lithosphere ( (cid:3) 120 km in the south) and mid-lithosphere sources to give a composition intermediate between magmas from the northern Lesser Caucasus and NW Iran. The mid-lithosphere magma source has a distinct composition compared with the base of the lithosphere, which is argued to be the result of the increased retention of metasomatic components in phases such as apatite and amphibole, which are stabilized by lower temperatures prior to magma generation.


INTRODUCTION
The Arabia-Eurasia collision zone is one of the very few places on Earth where it is possible to study active volcanism associated with a continent-continent collision event. The geodynamic processes that drive volcanism remain unresolved, with numerous and sometimes conflicting models for its origin (Pearce et al., 1990;Keskin, 2003;Gö gü s¸& Pysklywec, 2008;Faccenna et al., 2013). Most of the volcanism in the region is located on the Anatolian-Armenian-Iranian Plateau, a broad uplifted region to the north and east of the Arabian foreland, with an average elevation of over 2 km (S¸engö r et al., 2008;Priestley & McKenzie, 2013). The lithospheric structure of this plateau is considered to show a strong contrast between very thin mantle lithosphere in the west (below Eastern Anatolia), and very thick mantle lithosphere in the Zagros Core to the SE (Fig. 1), below western Iran (Priestley et al., 2012).
The Lesser Caucasus mountains sit close to the edge of the Zagros Core region, and therefore close to a transition from thick to thin lithosphere. The region thus provides an opportunity to look at the influence of lithospheric thickness on the geochemistry of collisionrelated magmas. Figure 1 shows the four volcanic highlands considered in this study, which form a NW-SE transect almost orthogonal to contours of lithospheric thickness from Priestley et al. (2012), which increase from NW to SE. It should be noted that the resolution on these lithospheric thickness estimates is limited by the 30-50 km vertical resolution of the seismic tomography data (McKenzie & Priestley, 2008). The Priestley et al. (2012) model is used because it is likely to give better local resolution in the Lesser Caucasus region than other global lithospheric thickness studies (Priestley & McKenzie, 2006, 2013. Volcanic rocks from the SE of this transect, where the lithosphere is thought to be thicker, are known to be more potassic than volcanic rocks from the NW (Meliksetian, 2013). This study provides the first complete geochemical dataset for volcanic rocks from the SE of the Lesser Caucasus. This dataset includes a complete range of compositions from basanite to rhyolite, which are used to evaluate the extent to which crustal contamination is an important component of magma petrogenesis. We then compare the geochemistry of the more primitive mafic volcanic rocks between the NW and SE of the Lesser Caucasus, and the mechanisms by which thicker lithosphere in the SE might influence the geochemistry of volcanic rocks found there. To indicate how thicker mantle lithosphere might influence the composition of magmas, volcanic rocks from a region of very thick lithosphere (>200 km) in NW Iran ( Fig. 1; Allen et al., 2013) are used as an end-member comparison of melting in a thick lithosphere regime.

Pre-collision geological history
The evolution of the Arabia-Eurasia collision zone, as a part of the Alpine-Himalayan orogenic belt, is a consequence of the closure of the Neotethys Ocean (Rolland, 2017, and references therein). The pre-Miocene geology of the territory of Armenia and the Lesser Caucasus (Fig. 2), like much of the interior of the Anatolian-Armenian-Iranian Plateau, is a complex amalgamation of a series of terranes (microplates), which accreted to the Eurasian continental margin during the closure of the Tethys Ocean (Hosseinpour et al., 2016;Rolland, 2017). The north and east of Armenia include rocks associated with the Mesozoic to early Cenozoic volcanic arc of the Lesser Caucasus (Mederer et al., 2013), an eastern continuation of the Pontides of Anatolia. The Pontide and Lesser Caucasus arcs together formed the active southern margin of the Eurasian Plate (Yilmaz et al., 2000). Figure 3 illustrates the various stages of  Neill et al. (2015). Red triangles denote the locations of major composite volcanoes. The locations of volcanic highlands that are used for the geochemical comparison that forms the focus of this study are shown by the coloured fields: 1, Shirak and Lori (Neill et al., , 2015; 2, Gegham (I. P. Savov, unpublished data); 3, Vardenis; 4, Syunik. New data for this study are from volcanic highlands 3 and 4. These volcanic highlands are close to parallel to major tectonic boundaries in the collision zone. Volcanic rocks from NW Iran referred to in the text are from the region labelled at the bottom of the map. Contours of lithospheric thickness from Priestley et al. (2012) are shown in red; numbers give lithospheric thickness in kilometres. The 'Zagros Core' refers to the region of maximum lithosphere thickness in the collision zone. closure of the Neotethys Ocean. During the Mesozoic there were probably several subduction zones that contributed to the closure of the northern and southern Neotethys basins (Fig. 3, panels 1 and 2; Galoyan et al., 2009;Rolland et al., 2010Rolland et al., , 2012Sosson et al., 2010;Karao glan et al., 2013Karao glan et al., , 2016Mederer et al., 2013;Topuz et al., 2013aTopuz et al., , 2013bTopuz et al., , 2014Hä ssig et al., 2015).
The numerous subduction zones illustrated in Fig. 3 would have added slab-derived material to large parts of the mantle wedge beneath the present-day collision zone. It has been suggested that this signature of a subduction-modified mantle has been widely inherited by the post-collisional volcanic rocks in the Arabia-Eurasia collision zone (Pearce et al., 1990;Keskin et al., 1998;Keskin, 2003;Allen et al., 2013;Oyan et al., 2017).

Post-Miocene volcanism and tectonics
The widespread volcanism occupies a wide zone across the collision zone ( Fig. 1), arguing against any single subducting slab driving magma generation. Elevated topography of 1-2 km is observed across the collision zone; however, the deep lithospheric structure that might isostatically support this elevated topography shows strong contrasts between the east and west. In the west, a low seismic velocity anomaly in the uppermost mantle below the East Anatolian Plateau (Al-Lazki, 2003;Angus et al., 2006;Zor, 2008;Gö k et al., 2011;Koulakov et al., 2012;Skolbeltsyn et al., 2014) has been used to suggest that the mantle lithosphere is very thin under Eastern Anatolia, such that the high topography is not isostatically supported. Based on the fairly uniform crustal thickness ($40 km; Angus et al., 2006) and the very thin lithospheric mantle, both the current structure of this western part of the Plateau and the associated volcanism have been attributed to a slab break-off event (Keskin, 2003(Keskin, , 2007S¸engö r et al., 2003, 2008. Given the later onset of eruptions, and the more alkaline nature of magmas in the south (Pearce et al., 1990;Keskin et al., 1998), it has been postulated that the slab that broke off and thereby fuelled recent volcanism was a north-dipping South Neotethys slab below the Bitlis Suture (Keskin, 2003).
In contrast, relatively fast seismic velocities in the upper mantle under western Iran suggest that a thick lithospheric root isostatically compensates the high topography in the southeastern part of the collision zone, termed the 'Zagros Core' (Priestley & McKenzie, 2006, 2013Priestley et al., 2012). There is no evidence in NW Iran of the coupled spatio-temporal variability in volcanism seen in Eastern Anatolia, such that instead of slab break-off, it has been postulated that the volcanism is driven by small-scale lithospheric dripping or amphibole breakdown melting Kaislaniemi et al., 2014;Lechmann et al., 2018).
As can be seen in Fig. 1, the Lesser Caucasus is located close to the boundary between these two tectonic regimes. The transect of volcanic highlands thus provides a key opportunity to understand the effect of lithospheric thickness on the composition of mantlederived magmas in a continental collision zone.

Geology and geochronology of collision-related volcanism in the Lesser Caucasus
Three major styles of volcanic activity can be observed in the Lesser Caucasus. The first of these is relatively low-volume eruptions from mostly small eruptive centres in monogenetic volcanic fields. In both the north and south of the Lesser Caucasus, this style of volcanism becomes increasingly dominant in the most recent eruptive products (Fig. 4;Connor et al., 2011;Neill et al., 2013). Second, large composite volcanoes are also found throughout the Lesser Caucasus. In the central Lesser Caucasus, such a volcano is Aragats (Connor et al., 2011), and in the south the smaller stratovolcanoes Tskhouk and Ishkhanasar (Fig. 4) dominate Fig. 3. Illustration of the progressive closure of the oceanic domains that separated Arabia from Eurasia during the Mesozoic, after Rolland (2017). SAB, South Armenian Block; SAB-BP, South Armenian Block-Bitlis Poturge terrane. The movement of the Bitlis Poturge terrane from Arabia to the South Armenian Block is not shown in these figures but occurs between panels 1 and 2. the landscape. Such large volcanoes (4090 m altitude, $70 km diameter in the case of Aragats), capable of generating many caldera collapse eruptions, are required to produce the ignimbrites that are widespread throughout the Lesser Caucasus. Finally, large fissure-fed 'flood basalt' style lava flows are found predominantly in some of the older ($3-2Á05 Ma) volcanic successions (Sheth et al., 2015).
Previous geochronological studies have constrained the ages of the volcanic deposits upon which this study focuses as being late Pliocene or younger. All of the samples analysed for geochemistry from the NW of the transect have ages of 2Á5 Ma or younger, based on correlations of the oldest 'valley series' lavas of Neill et al. (2015) with sediments containing mammalian fossils (Kharazyan, 1983), and a K-Ar age of 2Á5 Ma for one such 'valley series' lava in the Shirak province of NW Armenia in the northern Lesser Caucasus (Chernyshev et al., 2002). Volcanic rocks from further south were all collected from Quaternary volcanic highlands, which form monogenetic volcanic fields and so have an even more restricted age range. A limited number of ages from the Syunik volcanic highland (number 4 in Figs 1 and 2) suggest that volcanism ranges from 1Á3 to 0Á11 Ma, based on two Ar-Ar ages of pumice layers in diatomaceous sediments (Joannin et al., 2010), along with K-Ar ages of local lava flows that overlie the sediments . Archaeological evidence from 14 C dating of petroglyphs and burial places around the youngest lavas in Syunik suggests that volcanism may extend to within the last 5 kyr (Karakhanian et al., 2002). These ages are very similar to the estimated age span of the Gegham volcanic highland (highland 2 in Figs 1 and 2) of 1Á2-0Á02 Ma (Lebedev et al., 2013). In summary, the complete age span of the volcanic rocks studied for geochemistry in this study is <3 Myr, and it is likely that the present-day structure of the lithosphere can be used as an interpretive framework for all samples.

SAMPLING
Whereas recent work (Neill et al., , 2015 has provided important new geochemical data for collision-related volcanism in the northern Lesser Caucasus, our study presents the first comprehensive dataset for mafic volcanic rocks from the southern part of the Lesser Caucasus (Table 1). All samples are from the Vardenis and Syunik volcanic highlands (highlands 3 and 4 in Figs 1 and 2). Table 1 gives the coordinates for all the mafic samples in the Vardenis and Syunik volcanic highlands. Figure 4 shows the newly completed

Major and trace elements
Sample preparation was done at the University of Leeds. Rock samples (60-100 g) were crushed in a TEMA agate mill. The agate was cleaned extensively between the crushing of each sample, including precontamination of the agate by the crushing of 50 g of sample, which was then discarded. Bulk-rock major and trace element analysis of samples from this study was  Latitude:  40Á0905  40Á14205  40Á0768  40Á07712  40Á10283  40Á02458  40Á08497  40Á05927  40Á07163  Longitude:  45Á414  45Á6192  45Á51587  45Á5256  45Á39072  45Á83628  45Á83638  45Á79923  45Á78805  Elevation (m):  2257  2101  2340  2280  2188  2347  2375  2538 2539 Lab. analysis:* RH RH ACME RH ACME ACME ACME ACME ACME  done at ACME Labs by Bureau Veritas minerals, in Vancouver, Canada, and at Royal Holloway University, London. For analysis at ACME Labs, samples were heated to 1000 C to determine loss on ignition (LOI), and then fused in a platinum-gold crucible with a lithium tetraborate flux. The resulting fusion beads were then analysed by X-ray fluorescence (XRF) for major elements. Two internal ACME standards reproduced expected values to better than 3% (for all oxides >1 wt %). Trace element concentrations were determined by inductively coupled plasma mass spectrometry (ICP-MS) on the fused beads after digestion in nitric acid. Analysis of Ni concentrations involved digestion in aqua regia at 95 C. Two internal ACME standards reproduced expected values to better than 10% for trace elements. Values for internal and external standard data are shown in Supplementary Data A.
Analysis at Royal Holloway was by inductively coupled plasma optical emission spectroscopy (ICP-OES) for major elements and some high-abundance trace elements (Sr, Zr, Ni) and ICP-MS for low-abundance trace elements. Major element analyses followed the methods described by Walsh et al. (1981) and Garbe-Schö nberg (1993). The relative standard deviation on the external standards NIM-G, NIM-L, BHVO-1, RGM-1 and STM-1, as well as internal standards, was better than 5% (for all oxides >1 wt %). Trace elements were analysed after HNO 3 -HF-HClO 4 digestions. Prior to analysis, samples were spiked with 5 ng ml -1 of indium (In) and rhenium (Re) for internal standardization. For analytical quality control, five international   Nd/ 144 Nd are reported as the standard error, 2r (i.e. two standard deviations about the mean using 200-240 measurements). External precision (2r) for Sr and Nd isotopic ratios from successive replicate measurements of primary standards was better than 35 ppm for the NIST SRM-987 International Reference Standard ( 87 Sr/ 86 Sr ¼ 0Á710260 for nine runs averaged, with a standard deviation of 1Á1 Â 10 -5 ) and better than 25 ppm for the La Jolla Nd International Reference Standard ( 143 Nd/ 144 Nd ¼ 0Á511842 for 11 runs averaged, with a standard deviation of 2Á5 Â 10 -5 ). USGS standard BHVO-1 was also run as a validation material throughout the  Ba  1102  1176  1211  1152  1181  1180  1363  1093 1705 Sr  1585  2358  2576  1827  2150  2161  2198  2325  2812  Pb  15Á7  1 Zr  192  198  198  201  204  217  244  166
ubiquitous and normally occurs as glomerocrysts ( Fig. 5f). Some samples have abundant olivine microphenocrysts ( Fig. 5d) and phenocrysts (up to 600 lm), although olivine is generally restricted to the most mafic samples (low-silica trachybasaltic andesite or more mafic). Amphibole phenocrysts are common (up to 2 cm), and occur across the compositional range ( Fig. 5a). In several samples amphibole is absent, or partially resorbed (Fig. 5c), owing to disequilibrium conditions prior to eruption. The apatite phenocrysts can reach up to 2 mm in size in some cases (Fig. 5b). The common occurrence of hydrous minerals, such as amphibole and apatite, is noteworthy in comparison with the volcanic rocks from the northern and central Lesser Caucasus, where such mineral phases are less common, certainly in the most mafic samples (Connor et al., 2011;Neill et al., 2013). Typical groundmass in mafic volcanic rocks from Syunik is composed of plagioclase, clinopyroxene, and Fe-Ti oxides, with volcanic glass in some scoria samples. It should be noted that plagioclase is present only in the groundmass and is not a phenocryst phase. More evolved high-silica trachybasaltic andesites, trachyandesites and trachytes from Syunik show the appearance of abundant plagioclase phenocrysts. In some cases, these plagioclase crystals (up to 2 mm) show evidence of multiple stages of crystallization with distinct core to rim zonation (Fig. 5e). Amphibole also becomes a more common phenocryst phase, whereas clinopyroxene is less dominant. The groundmass is increasingly dominated by plagioclase. Rhyolites are common in the northern part of Syunik, where they commonly form obsidian flows.
In Vardenis all samples are trachybasaltic andesites or more evolved, and plagioclase is always a phenocryst phase. In one trachybasaltic andesite (Table 2, sample 6.3.15), biotite is also a phenocryst phase. Rhyolites from Vardenis commonly have biotite, plagioclase, quartz and potassium feldspar as phenocrysts, and a groundmass composed of quartz and potassium feldspar.
From the occurrence of different minerals in the range of rock types we sampled across the Syunik and Vardenis volcanic highlands (34 samples; Table 2), it is possible to suggest a general order of crystallization for both areas, which would probably have been: olivine þ Fe-Ti oxides 6 apatite; clinopyroxene þ Fe-Ti oxides 6 olivine 6 amphibole; clinopyroxene þ plagioclase þ Fe-Ti oxides 6 amphibole 6 phlogopite; plagioclase þ potassium feldspar þ Fe-Ti oxides 6 phlogopite. The predominant focus of this study is on the mechanism of mantle melting, which requires the effect of fractional crystallization to be minimized. Only primitive samples with <54 wt % SiO 2 and >4 wt % MgO are used to investigate questions of magma petrogenesis. These samples will probably all have only fractionated mafic minerals, with no significant feldspar fractionation.

Major element characteristics
Across the Lesser Caucasus there is a great diversity in the compositions of the volcanic rocks within each volcanic highland, with most showing a complete compositional range from basalt to rhyolite (Fig. 6). Rhyolites are also present in the northern Lesser Caucasus (Karapetian et al., 2001), although they were not sampled for the Neill et al. (2013)   extreme, have more alkaline compositions compared with those from the north, with the largest Na 2 O and K 2 O variations between north and south in the most mafic samples (Fig. 6). Southern Lesser Caucasus samples are also more potassic when compared with their northern counterparts (shoshonitic versus calc-alkaline series; Fig. 7a). As well as being more alkaline, mafic southern Lesser Caucasus samples extend to lower SiO 2 contents.

(d) c) c) c) c) c) c) c) c) c) c) c) c) ) ) c) ) ) ) c) c) (c)
On some major element variation diagrams, such as MgO versus SiO 2 (Fig. 7b), the trends are fairly similar for northern and southern Lesser Caucasus samples. However, other elements, notably P, display a significant geographical gradient in concentration for the most mafic samples, from 0Á4 wt % P 2 O 5 in the north to as high as 1Á6 wt % in the south (Fig. 7c).

Trace element characteristics
Mid-ocean ridge basalt (MORB)-normalized trace element patterns all show profiles typical of subductionrelated volcanic rocks, with positive anomalies for Ba, K, Pb and Sr and negative anomalies for the high field strength elements (HFSE) Nb, Ta and Ti (Fig. 8). Superimposed on this is an additional enrichment in incompatible trace elements, in particular the light rare earth elements (LREE; Fig. 9), Sr, Ba and P, which becomes increasingly pronounced to the south of the transect. Owing to the LREE enrichment, as well as a more moderate depletion in heavy rare earth elements (HREE), REE profiles become increasingly steep towards the south (Fig. 9), with CI-normalized La/Yb ratios of five in the northern Lesser Caucasus compared with values as high as 40 in the south.
These additional enrichments do not resemble those of typical intra-plate ocean island basalts (OIB), which show uniform enrichments in all the most incompatible trace elements, rather than the larger enrichments seen in the LREE relative to the HFSE (Fig. 8e), as exemplified by Ta and La in Fig. 7e and f, respectively. These additional enrichments are seen most clearly in volcanic rocks from NW Iran (Fig. 8e), a region of very thick lithosphere-up to 220 km (Priestley et al., 2012). In Fig. 8a-d we show that the variations in the composition of basalts within individual volcanic highlands are small, relative to the variations in composition between the highlands. Both the ubiquitous subduction ('arc') signatures and additional enrichment are highlighted in Fig. 10. All samples plot above the mantle array, typical for rocks from volcanic arcs. However, the samples from the southern Lesser Caucasus have higher Th/Yb and Ta/Yb relative to their northern counterparts (Fig. 10). Volcanic rocks from the thick lithosphere Zagros Core region of NW Iran again plot as a compositional end-member.
Given that amphiboles are relatively common in the southern Lesser Caucasus rocks, it should be noted that low Y (amphibole fractionation) trends dominate the Y versus SiO 2 plot for these rocks (Fig. 7d). Samples from further north, in which anhydrous mineral assemblages are more common, show both high Y (anhydrous assemblage) and low Y (hydrous) fractionation trends.

Sr-Nd isotope systematics
The Sr-Nd isotope compositions of the southern Lesser Caucasus samples are shown in Table 3. 87 Sr/ 86 Sr varies from 0Á7043 to 0Á7047, whereas eNd varies from þ2Á2 to þ4Á2. All mafic volcanic rocks from the Lesser Caucasus plot along the mantle array, with compositions more enriched than normal (N)-MORB but more depleted than Bulk Earth (Fig. 11). In general, those samples from the SE have higher (more radiogenic) 87 Sr/ 86 Sr, and less radiogenic eNd than those samples from the NW, defining a steep gradient on the Sr-Nd isotope diagram (Fig. 11). A few samples plot away from this trend to higher 87 Sr/ 86 Sr (see Crustal contamination section of Discussion). Lesser Caucasus volcanic rocks are isotopically more depleted than volcanic rocks formed above very thick lithosphere in NW Iran (>200 km; Priestley et al., 2012). They are also more depleted than volcanic rocks from the Lake Van area, which commonly display significantly more variable Sr-Nd isotope compositions (not shown), probably because those magmas interacted extensively with continental crust during ascent (Pearce et al., 1990;S¸en et

Crustal contamination
Volcanic rocks produced within the thicker lithosphere of the southern Lesser Caucasus have higher concentrations of incompatible trace elements, more radiogenic 87 Sr/ 86 Sr, and lower 143 Nd/ 144 Nd ratios. All of these features could be produced by crustal contamination during magmatic evolution, as has been suggested for various parts of Eastern Anatolia, where assimilation of radiogenic ancient continental crust gives 87 Sr/ 86 Sr ratios of up to 0Á7065 and marked variation of these ratios with SiO 2 (Pearce et al., 1990;Keskin et al., 2006). All Sr-Nd isotope ratios for the Lesser Caucasus volcanic rocks show small variations compared with what would be expected if the magmas had been contaminated by continental crust. The crust in the Lesser Caucasus is composed of a mixture of felsic metamorphic basement (South Armenian Block), arc volcanic rocks and Mesozoic-Cenozoic sediments (Fig. 2). Only the arc-related volcanic rocks would have similar isotope compositions to the collision-related magmas, such that assimilation of any of the other lithologies would alter, often dramatically, the isotope composition of the host magma. The majority of samples show insignificant variability in 87 Sr/ 86 Sr and 143 Nd/ 144 Nd with SiO 2 within individual volcanic highlands (Fig. 12), suggesting that isotope ratios are not being modified significantly during magma evolution and storage in the crust prior to eruption. The lack of isotopic variability with SiO 2 content suggests that assimilation of South Armenian Block crust (Fig. 2), which with a basement 87 Sr/ 86 Sr ratio of $0Á7303 (Baghdasaryan & Ghukasyan, 1985) should be easily identified, is unlikely. Similarly, assimilation of sedimentary material is also unlikely, given typical Tethyan (Mesozoic) flysch 87 Sr/ 86 Sr of 0Á7112 (Prelevi c et al., 2008).
One sample from the Syunik volcanic highland (1-4A-08), which was sampled from a scoria cone containing large (centimetre-sized) felsic xenoliths entrained within the trachybasaltic andesite scoria, has trace element and isotopic characteristics indistinguishable from the other mafic samples. The Sr-Nd isotope composition of one of the xenoliths (sample 1-4B-08; Table 3) is also shown in Fig. 12, and is only slightly elevated in 87 Sr/ 86 Sr (0Á7049) above that of the Quaternary basalts, suggesting that an unrealistically high degree of assimilation would have to occur for the composition of the magma to be significantly affected. The 143 Nd/ 144 Nd (0Á5128) is indistinguishable from Syunik basalts. It is likely that the felsic xenolith is cogenetic with the trachybasaltic andesite scoria host. It provides evidence that basaltic magmas in the southern Lesser Caucasus are interacting with felsic igneous rocks in the crust, which have a similar origin, rather than interacting with ancient crust, which would have strongly affected the Sr-Nd isotope ratios.
The curved patterns for major element variations versus SiO 2 , including MgO (Fig. 7b), P 2 O 5 (Fig. 7c) and Al 2 O 3 (not shown), suggest that the rhyolites are derived from extreme degrees of magmatic differentiation, rather than being crustal melts. These rhyolites also do not appear to be mixtures between primitive magmas and continental crust (Rudnick & Gao, 2003). The Rb and Sr contents of a rhyolite from Syunik are  Fig. 1, respectively. The selected samples are the three samples with the highest MgO content (wt %) in each volcanic highland. Data sources are as in Fig. 6. For comparison, the total range of 'basalts' (<52% SiO 2 ) from the Syunik volcanic highland is shown in (a) and (b) as the pale green field. The total range of Shirak and Lori basalts is shown as a gold field in (c) and (d). (d) also shows the composition of a Syunik rhyolite (diamonds) and bulk continental crust (stars) for Rb, K and Sr, to demonstrate whether these rhyolites are likely to be formed by assimilation and/or melting of continental crust. Composition of the continental crust is from Rudnick & Gao (2003). The total compositional variability is illustrated in (e), where averages from the geographical extremes (Shirak in the north and Syunik in the south) are compared. Also shown here is a sample from Kurdistan, NW Iran from Allen et al. (2013), formed from melting in a region of very thick lithosphere. Normalization factors are from Sun & McDonough (1989). The average composition of continental arc basalt is from Kelemen et al. (2003). actually more extreme relative to the basalts than average continental crust (Rb more strongly enriched, Sr more strongly depleted in the rhyolite; Fig. 8d). K 2 O is more enriched in the rhyolite than in the basalts, whereas continental crust has lower K 2 O than continental crust (Fig. 8d). The lack of evidence for extensive assimilation of continental crust in the petrogenesis of the rhyolites suggests that these processes are likely to be unimportant in the petrogenesis of the more primitive magmas.
It is also possible that the magmas may have been contaminated by Mesozoic-Paleogene arc crust (Fig. 2) on their ascent, as has been suggested for the northern Lesser Caucasus in some isolated cases (Neill et al., 2015). The similar Sr-Nd isotope compositions of these arc rocks (Mederer et al., 2013) to the collision-related magmas mean that assimilation could be 'cryptic', without obvious modification of Sr-Nd isotope ratios. Given the arc-like incompatible trace element geochemistry of Lesser Caucasus volcanic rocks, arc-crust assimilation is unlikely to significantly alter trace element compositions, although it is possible that such assimilation could explain some of the spread in the southern Lesser Caucasus data. The average composition of Mesozoic arc rocks from the southern Lesser Caucasus (Kapan Zone) is shown in Fig. 10 (Mederer et al., 2013). The Th/ Yb ratio of the arc rocks is elevated above the mantle array similar to the post-collisional Lesser Caucasus samples. However, the low Ta/Yb ratios of these arc rocks make them unlikely candidates to explain the more enriched composition of southern Lesser Caucasus magmas.
There are a few examples of elevated 87 Sr/ 86 Sr ratios in evolved samples above the background range. In most cases, it seems that the Nd isotopes are unaffected. Some ancient arc samples from the Kapan zone, southern Lesser Caucasus, have very high initial 87 Sr/ 86 Sr, but 143 Nd/ 144 Nd does not vary significantly in the same samples (Mederer et al., 2013). Extensive assimilation of such material would be required to explain the elevated 87 Sr/ 86 Sr of these samples. Worthy of note, however, are the very high Rb/Sr ratios (7-23) in these rhyolites-sometimes 1000Â greater than those of typical basalts (Table 1). Whereas decay of 87 Rb over the relatively short timeframe (<2 Myr) since cooling and crystallization of lavas will have negligible effects on the 87 Sr/ 86 Sr of basaltic rocks, such high Rb/Sr ratios mean that post-crystallization decay will have a significant effect for rhyolites (Fig. 12a). Indeed, minimum ages for the rhyolites of between 0Á6 and 1Á4 Ma are sufficient to give initial 87 Sr/ 86 Sr within the range observed for basalts. As such, even in these rhyolites, crustal contamination may have played an unimportant role in magma evolution. It should be noted that the lower Rb/ Sr ratio of the felsic xenolith (0Á12) means that its elevated 87 Sr/ 86 Sr ratio was probably present during emplacement of the host Quaternary mafic magma.
In summary, assimilation of crustal rocks with markedly distinct geochemistry is unlikely, and although assimilation of igneous rocks with similar petrogenetic origins cannot be precluded, such processes would be incapable of driving the enrichment in incompatible trace element concentrations and Sr-Nd isotope ratios observed from north to south in the volcanic rocks of the Lesser Caucasus. The lack of evidence for crustal contamination, along with the gradient in Sr-Nd isotope ratios in Fig. 11, requires that there must be variation in the composition of the mantle source.

Subduction modification of the mantle source in the Lesser Caucasus
All mafic samples show distinct negative Nb-Ta and Ti anomalies (Fig. 8), and positive spikes in Ba, K, Pb and Sr in normalized trace element patterns, such that the overall patterns, if not the absolute concentrations, are typical of arc magmas. In the absence of evidence for widespread crustal contamination, it is likely that these features are inherited from the mantle source. Collision between Arabia and Eurasia was preceded by subduction of various Tethyan ocean basins along several convergent margins (Fig. 3). These subducting slabs would probably have contributed slab material to the mantle below the collision zone.
Following the approaches of Hofmann (2003) and Turner et al. (2017), Fig. 13 illustrates how the ratios Ba/ La, Ce/Pb, Sr/Nd and Nb/U in mafic samples (<54 wt % SiO 2 , >4 wt % MgO) vary with latitude. These ratios show little variation within all MORB and OIB, but are much more variable in volcanic arc rocks owing to the subducted slab contribution, such that they can be used here as proxies. In all of the plots in Fig. 13 these ratios deviate from those of MORB or OIB, demonstrating the presence of a slab component in their mantle source. These ratios do not show consistent variability between the volcanic highlands ( Fig. 13), although Sr/Nd ratios are slightly elevated to the south (Fig. 13c) and the Nb/U ratio reaches a minimum (highest slab contribution) in the central Lesser Caucasus. The lack of a consistent trend in these ratios between volcanic highlands means that the slab contribution is likely to be fairly uniform.
Although mixing between a depleted mantle source and the likely composition of subducted sediment (e.g. Tethyan flysch; Prelevi c et al., 2008) comes close to explaining the Sr and Nd isotopic composition of the least enriched samples in Fig. 11, such mixing is unable to explain the trend towards the higher 87 Sr/ 86 Sr and lower 143 Nd/ 144 Nd seen in the SE. This enrichment in the southern Lesser Caucasus must therefore be driven by some other process. Pre-collision subduction events have probably imparted a subduction signature on the mantle source across the Lesser Caucasus. However, the geochemical gradient between volcanic rocks of the northern and southern Lesser Caucasus cannot be explained by differences in the composition or size of the slab contribution to the mantle source.

Modelling the conditions of mantle melting
It is possible to investigate how a thicker lithosphere affects the conditions of melting by using the approach of Shaw (2005) to forward model the composition of the samples using a non-modal batch melting model (Table 4), with the equation  Pearce (1983). Lesser Caucasus data sources are as in Fig. 6. NW Iran data are from Allen et al. (2013). Fractional crystallization (FC) vector: from basalt to andesite uses a 50% amphibole þ 50% clinopyroxene assemblage and basaltic melt partition coefficients; from andesite to rhyolite uses a 50% plagioclase þ 50% amphibole assemblage and andesite melt partition coefficients. The vector has a starting composition as the most mafic sample. Partition coefficients for Ta are: clinopyroxene-basalt ¼ 0Á017; amphibole-basalt ¼ 0Á05; amphibole-andesite ¼ 0Á21; plagioclase-andesite ¼ 0Á03. For Th: clinopyroxene-basalt ¼ 0Á007; amphibole-basalt ¼ 0Á05; amphibole-andesite ¼ 0Á16; plagioclase-andesite ¼ 0Á01. For Yb: clinopyroxene-basalt ¼ 0Á28; amphibole-basalt ¼ 0Á59; amphibole-andesite ¼ 1Á25; plagioclase-andesite ¼ 1Á25. All partition coefficients are from the GERM database (earthref.org/GERM/). The average compositions of Mesozoic arc rocks from the Kapan zone arc rocks are from Mederer et al. (2013). The grey box labelled BCC is the composition of bulk continental crust from Rudnick & Gao (2003).
where C l is the concentration of an element in the liquid, C 0 is the concentration of that element in the source, D is the element's bulk partition coefficient, F is the fraction of melting, and P represents the partitioning of the element into the melt according to the proportion in which each mineral enters the melt. None of the samples in the region are in equilibrium with mantle olivine, so all must have undergone some fractional crystallization between last equilibration with the mantle and eruption. This is accounted for by assuming 8% fractional crystallization of an assemblage composed of 90% olivine and 10% spinel [following the approach of Shaw (2005)]. From equation (1) we can see that there are three parameters, each of which can be varied to generate the trace element concentration of the most primitive southern Lesser Caucasus samples. First, the fraction of melting (F) could vary, affecting the concentration of all incompatible trace elements. Second, if the modal mineralogy is changed then the partitioning of elements between source and melt [D and P in equation (1)] will change. Third, the concentration of elements in the source rock (C 0 ) could be changed. Our approach is to take the equivalent model of Neill et al. (2015) for the most geochemically depleted samples from the northern Lesser Caucasus and iteratively vary two of these parameters to attempt to reproduce the average composition of primitive basalts in the south. Variations in the mineralogy of the melt source are simplified to spinel versus garnet peridotite melting. Spinel peridotite modal mineralogy and melting proportions are from Shaw (2005) (Jacobsen & Wasserburg, 1980). *Samples from the east side of the Sevan-Akera suture (see text for discussion).
respectively, whereas for the garnet peridotite they are from Thirlwall et al. (1994) and Allen et al. (2013), respectively (Table 4). Varying the fraction of melting is self-explanatory. Remaining disparities between the model and observed trace element concentrations are likely to be due to differences in the composition of the mantle source. The different stages of melt modelling are shown in Table 4. REE chemistry is used to obtain a qualitative understanding of the changing conditions of melting between the northern and southern Lesser Caucasus (Fig. 14). In Fig. 14a, two vectors are plotted, which show how the composition of the magmas should change with the degree of melting versus source mineralogy (presence or absence of garnet) and/or changes in source composition. Davidson et al. (2013) showed that REE partitioning between clinopyroxene and basaltic magma means that melting curves for peridotite will be very steep on a Dy/Dy* versus Dy/Yb plot. The assumption that clinopyroxene is a residual phase is reasonable given the low degrees of melting (3% or less) previously estimated for the northern Lesser Caucasus (Neill et al., 2015). Amphibole is also likely to be an important residual phase given the positive correlation between Dy/ Dy* and Ti/Ti* (the size of the Ti anomaly in Fig. 8; not shown here). It is unclear from Fig. 14a whether garnet in the mantle source, or a LREE-enriched source is responsible for the high Dy/Yb ratio in southern Lesser Caucasus samples. Both should give vectors close to horizontal in Fig. 14a, as both changes will steepen the REE profile rather than changing its curvature (Davidson et al., 2013). The southern Lesser Caucasus samples have lower HREE abundances than samples from the north (Fig. 9). HREE all behave incompatibly during melting of spinel peridotite (Table 4), such that lower degrees of melting as suggested by the lower Dy/ Dy* ratio should only increase the concentration of HREE. Therefore, garnet as a residual phase during melting must explain at least some of the increase in Dy/Yb in Fig. 14a. In Fig. 14b, the northern Lesser Caucasus samples sit close to the spinel peridotite melting curve, but the southern Lesser Caucasus samples are intermediate between the spinel and garnet melting curves, suggesting melting of a mixed source involving both garnet and spinel peridotite. This could be a consequence of polybaric melting across the spinel-garnet Fig. 11. eNd vs 87 Sr/ 86 Sr for volcanic rocks from the Lesser Caucasus. Bulk Earth value and mantle array are after Rollinson (1993). NW Iran data are from Allen et al. (2013); this region has very thick lithosphere, estimated to be >200 km by Priestley et al. (2012)  transition at $75 km depth (Robinson & Wood, 1998), suggesting melting at greater depths in the south.
Returning to the partial melting model, the melt fraction was estimated on the basis that the Nb concentration (Table 4) is not affected by the mineralogy of the source rock, and assuming that its concentration in the source rock is constant along the transect. Utilizing this assumption gives 3% melting in the north, versus 1% melting in the south (Fig. 15). The proportions of garnet peridotite and spinel peridotite that contribute to the total mantle source are estimated from the Hf and Yb concentrations. Both elements are assumed to have constant concentrations in the mantle source because HFSE and HREE would be less affected by metasomatic events that could alter the source composition. The best fit to the Hf and Yb data is a magma source that is 65% garnet peridotite and 35% spinel peridotite (Fig. 15).
Across the Lesser Caucasus, as lithosphere thickness increases, the degree of melting decreases and melting occurs at greater depths, with garnet as a residual phase for a significant portion of the melting interval. However, the model is still unable to explain the high concentrations of several of the LREE (La, Gd and Dy; Fig. 15, black line), and requires changes in the composition of the mantle source with lithospheric thickness as well.
To understand how melting is occurring, it is important to understand where magma is forming with respect to the lithospheric structure. For this, we have used major element thermobarometry to estimate the conditions of last equilibration between magmas and the mantle.

Pressure and temperature of melting
To calculate these intensive parameters, we use the parameterizations of Plank & Forsyth (2016), based on the major element chemistry of primitive magmas, building on the work of Lee et al. (2009). These calculations are for magmas produced from a peridotite mantle source at pressures below 3 GPa (Plank & Forsyth, 2016). The temperature (T) dependence of Mg partitioning between olivine and melt (Roeder & Emslie, 1970) and the pressure (P) dependence of silica activity in melts co-saturated in olivine and orthopyroxene (Carmichael et al., 1970) are exploited to give the equations (2) All major element oxides are calculated in mol % as described by Lee et al. (2009), except that they are calculated on an anhydrous basis. The two terms DT H2O and DT CO2 DT CO2 account for the lower temperature and higher pressure melting in the presence of volatiles as follows: DT H2O ¼ 40Á4ðH 2 OÞ À 2Á97ðH 2 OÞ 2 þ 0Á0761ðH 2 OÞ 3 (4) where H 2 O and SiO 2 are in wt %. To use these equations, the composition of the primary melt must be known. Given that no samples from the Lesser Caucasus have Mg# > 67, it is clear that all the studied samples have undergone some degree of fractional crystallization. If only olivine has fractionated, then it is a simple process to add olivine incrementally to the melt until it is in equilibrium with mantle olivine of Fo 90 (Lee et al., 2009). If, however, other phases, such as amphibole or clinopyroxene, have also fractionated, then this calculation becomes non-trivial. In a study on the collision-related volcanism of Anatolia, McNab et al. (2018) took 8Á5 wt % bulk-rock MgO as a reasonable lower limit for magmas that have fractionated only olivine. Only two such samples exist in the Lesser Caucasus, and it is likely that they have both been affected by olivine accumulation rather than being truly primitive melts. Instead, the complete sample suite from the Shirak and Lori, and Syunik volcanic highlands (1 and 4 in Fig. 1) is used to project back to the likely composition at 8Á5 wt % MgO of two end-member primitive magmas for the Lesser Caucasus. This correction is shown in Supplementary Data B.
To calculate the pressures and temperatures of melting, two magma composition parameters must be constrained: Fe 3þ /RFe and water content. The water content of the primary magma has a large effect on the temperature estimate through equation (4) ($25 C per wt %), but a smaller effect on the pressure estimate ($30 MPa per wt %). The presence of hydrous mineral Melting models La Gd Dy Yb Zr Hf Nb 1. 3% melting spinel peridotite (Neill et al., 2015) Initial The source modes for spinel and garnet peridotite are from Shaw (2005) and Thirlwall et al. (1994) respectively. Melt modes are from Allen et al. (2013) and Neill et al. (2015). The starting composition of the mantle source is the same as that used by Neill et al. (2015). Partition coefficients are from Ionov et al. (2002), except for those for garnet, which come from McKenzie & O'Nions (1991) and Adam & Green (2006). The most primitive Syunik sample used for comparison is 8-5-15. OL, olivine; OPX, orthopyroxene; CPX, clinopyroxene; AMPH, amphibole; SP, spinel; GRNT, garnet.
phases in mafic rocks (as opposed to only in more evolved samples), only in the southern Lesser Caucasus (Fig. 5), suggests higher water contents in the south. Based on previous studies of amphibole peridotite xenoliths (Thirlwall et al., 1994;Ionov & Hofmann, 1995), it is assumed that the mantle source contains between 2 and 10 modal % amphibole, which is assumed to contain 2 wt % H 2 O, such that the mantle source could have up to 0Á2 wt % H 2 O. Using Ce partition coefficients for water during mantle melting (Thirlwall et al., 1994;Dixon et al., 2002;Ionov et al., 2002), and the melting models shown in Table 4, gives 2-7 wt % H 2 O for southern Lesser Caucasus primary magmas, and 1-4Á6 wt % for the north.
As only ferrous iron substitutes in olivine, a higher Fe 3þ /RFe increases the apparent melt Mg#, and so reduces the amount of olivine addition required to produce the primary magma, the MgO content of the primary magma, and thus the calculated pressure and temperature. The Fe 3þ /RFe ratio varies from a minimum of 0Á1 in MORB (Cottrell & Kelley, 2011) to as high as 0Á3 in some arc basalts (Brounce et al., 2014). Given the ubiquitous arc-type geochemical signatures in Lesser Caucasus volcanic rocks, we take 0Á25 as a conservative estimate of the Fe 3þ /RFe ratio (Brounce et al., 2014).
It should be noted that equations (2) and (3) will give meaningful estimates of the pressure and temperature of melting only if the mantle source is peridotitic, because the parameterizations assume the melt is saturated in both orthopyroxene and olivine (Plank & Forsyth, 2016). The use of these equations on magmas derived from a pyroxenite source would yield meaningless results. Figure 16 shows that for Lesser Caucasus and NW Iran samples, pyroxenite is not a major component of the mantle source. Pyroxenite partial melts have much higher Ni/MgO ratios. As magmas fractionate, those derived from a peridotite source should evolve along a trajectory below the dashed line in Fig. 16, whereas those from a pyroxenite source would plot above the line.  Davidson et al. (2013). The two arrows depict the expected vectors from changing the fraction of melting or source mineralogy or composition. The vector for a lower fraction of melting is after Davidson et al. (2013, fig. 4), based on melting of an olivine-pyroxene-amphibole-dominated mantle source. The presence of garnet in the magma source should follow a horizontal vector after Davidson et al. (2013, fig. 5). (b) Dy/Yb vs La/Yb, melting curves based on the modal mineralogies, melting modes, and partition coefficients given in Table 5. Fig. 15. N-MORB-normalized concentrations of Nb, La, Zr, Hf, Gd, Dy and Yb for the average of mafic samples from Syunik volcanic highland (green line) and several non-modal batch melting models that attempt to explain these trace element concentrations. Melting models: spinel peridotite at 3% melting (yellow dashed line), 1% melting of spinel peridotite (dashed black line), 1% melting of garnet peridotite (black continuous line) and 1% melting of a source composed of 65% garnet peridotite and 35% spinel peridotite (bold dash-dot line). The blue dashed line represents a magma derived from the same 65% garnet peridotite source, with 4% apatite dissolved in the melt. Apatite composition is for magmatic apatite from Western Turkey of Prelevi c et al. Based on our water content estimates, southern Lesser Caucasus magmas last equilibrated with the mantle at between 1198 and 1292 C, whereas northern Lesser Caucasus magmas last equilibrated at between 1174 and 1267 C. Given the observation that amphibole is likely to be a residual phase in the mantle source (see previous section), it is likely that melting was occurring at temperatures close to the amphibole dehydration solidus, because amphibole is likely to completely break down within a few tens of degrees centigrade of crossing the dehydration solidus (Green & Falloon, 2005;Mandler & Grove, 2016). As such, it is likely that the minimum temperatures of these ranges are the more realistic. Using these temperatures gives pressure estimates of 2Á1 GPa ($75 km) in the southern Lesser Caucasus, and 1Á2 GPa ($45 km) in the north. These estimates are shown in Fig. 17, along with the position of the amphibole dehydration solidus. Southern Lesser Caucasus magmas are produced deeper, but at similar temperatures to those in the north.
The parameterizations first developed by Lee et al. (2009) were designed to expand the applicability of basalt melt geothermobarometry beyond mid-ocean ridge systems to any setting involving the melting of terrestrial peridotite. Both our modelled major element compositions and the pressure and temperature calculated fall within the experimental range of the dataset of Lee et al. (2009), and to this extent our approach is justified.
However, the presence of metasomatic phases such as amphibole in the mantle source means that the pressures and temperatures calculated here must be interpreted with caution. Alkaline volcanic rocks from monogenetic volcanic fields in Western Mexico, geochemically similar to those from volcanic highlands in the Lesser Caucasus, have very high Fe 3þ /RFe ratios and whole-rock Mg# in excess of the $72 value normally assumed to be in equilibrium with mantle peridotite (Carmichael et al., 1996). The effects of a potential underestimation of both Fe 3þ /RFe and the amount of olivine addition required to produce the primary magma will tend to counteract each other, where the former would give an overestimation of temperature, and the latter an underestimate. In the case of Western Mexico, both of these observations are ascribed to the incongruent melting of phlogopite (Carmichael et al., 1996), which does not seem to be a major metasomatic phase in the Lesser Caucasus mantle (see the section 'How does thicker mantle lithosphere influence the composition of the mantle source?' below).
Even if these issues are relevant to the Lesser Caucasus, the constraint provided by the position of the amphibole dehydration solidus shows the temperature estimates provided here to be reasonable, although some of the $50 C disparity with the dehydration solidus could be explained by these issues. The sensitivities of the pressure estimate to Fe 3þ /RFe (-0Á2 GPa per 0Á1 increase), primary melt Mg# (þ0Á5 GPa per 10% increase) and estimated temperature (þ0Á1 GPa per 100 C) are all small. The pressure estimate is primarily based on silica activity in a system co-saturated in olivine and orthopyroxene. Given the indications from Fig. 16, and reasonably successful trace element modelling of partial melting of a peridotite mantle source, it seems likely that these pressure estimates are robust.
Further constraints on the thickness of the crust are required to interpret these results. Crustal thickness is estimated using the formulations of Hu et al. (2017), which link the Sr/Y and (La/Yb) N ratios of intermediate magmatic rocks from continental collision zones with the crustal thickness. The basis of this technique is the Fig. 16. Ni/MgO vs SiO 2 in primitive Lesser Caucasus and NW Iran bulk-rock samples (<54% SiO 2 , >4% MgO), after Allen et al. (2013). Pyroxenite melting should generate compositions above the bold red dashed line (Sobolev et al., 2005). Primary melts of peridotite and pyroxenite (from Hawaii), after Sobolev et al. (2005). FC, fractional crystallization. Syunik, Shirak and Lori and NW Iran data are shown with the same symbols as in other figures. Fig. 17. Depth vs temperature of melting for the northern Lesser Caucasus and southern Lesser Caucasus. Also shown are the anhydrous peridotite solidus (1100 C þ 3Á5 C km -1 ) after Plank & Forsyth (2016), the wet solidus (amphibole present) after Green & Falloon (2005), and samples from East Anatolia with low K/Nb (red circles) from McNab et al. (2018). The pressures and temperatures of melting of NW Iran magmas were calculated for this study assuming Fe 3þ /RFe of 0Á25 and melt water contents of 7 wt %, as were used for the southern Lesser Caucasus estimate. (See text for discussion.) polarizing effects different fractionating mineral assemblages have on the Sr/Y and (La/Yb) N ratios in shallow versus deep storage reservoirs. To the nearest 5 km, we estimate crustal thickness as 45 km in Shirak and Lori, 55 km in Gegham, and 60 km in both Vardenis and Syunik. The depth of melting in the northern Lesser Caucasus is very similar to the $45 km Moho depth, possibly suggesting the presence of a very thin mantle lithosphere. It is instructive to compare this result with recent estimates of melting conditions in neighbouring Eastern Anatolia from McNab et al. (2018). East Anatolian magmas can be split into high K/Nb (>500) and low K/Nb types. Low K/Nb magmas generally have OIB-like geochemistry, and are considered to be directly derived from melting of the convecting mantle, which is apparently anomalously hot in the region, plotting to the right of the ambient mantle adiabat (Fig. 17), at a much higher temperature than any temperatures modelled for the northern Lesser Caucasus. This can be interpreted as northern Lesser Caucasus magmas being derived from the lithosphere, and not the convecting mantle.
Deeper melting in the south is associated with thicker crust, but magma generation is occurring at much shallower depths than the >100 km lithosphereasthenosphere boundary estimated by Priestley et al. (2012). The depth of melt equilibration is much deeper than the anhydrous solidus, confirming the need for a volatile-enriched lithospheric mantle root.
These P-T conditions of melting can be compared with estimates made for this study of the melting conditions of samples from NW Iran, based on the data of Allen et al. (2013). The P-T conditions of melting of NW Iran volcanic rocks were estimated using the same Fe 3þ /RFe and water contents as for southern Lesser Caucasus magmas (0Á25 wt % and 7 wt %, respectively). This is on the basis that the geochemistry of NW Iran magmas is similar to that of southern Lesser Caucasus magmas (Fig. 8), and both formed in a thicker lithosphere regime. Most of the magmas in NW Iran formed under similar conditions to southern Lesser Caucasus magmas, with a small subsidiary group of samples that formed at shallower depths, similar to northern Lesser Caucasus magmas, possibly representing magmas that re-equilibrated with the mantle at the base of the crust during their ascent to the surface.
The melting depth in both the southern Lesser Caucasus and NW Iran is significantly shallower than the estimated lithospheric thicknesses. One possible explanation for what these melting conditions represent could be a thermal maximum in a back-bent non-linear geothermal gradient (Fig. 17). In recently thickened orogenic lithosphere, this type of kinked profile is probably more realistic than the linear geothermal gradient more typical of a cratonic region (Mather et al., 2011).
It is worth noting that these differences in melting conditions probably reflect variations in the location of magma generation today, and are not the product of temporal changes in lithosphere structure. This is because, in Fig. 8, the trace element compositions of basalts from each volcanic highland are distinct, with each having a fairly narrow range of trace element compositions, despite all the volcanic highlands spanning an age range >1 Myr, and with the Gegham and Syunik highlands thought to have very similar timespans of activity (Joannin et al., 2010;Lebedev et al., 2013).

Geodynamic implications
Using constraints from geochemistry and the estimated pressures and temperatures of melting, it is possible to develop a model of how melting occurs under the Lesser Caucasus region. It seems clear that a majority of melting is taking place in the lithosphere for two reasons. First, all samples display an arc-type geochemistry. Subduction ceased around 35 Ma (Rolland et al., 2012), such that these arc signatures will probably be preserved only in the lithosphere. It is possible that the mantle lithosphere retains phases such as amphibole and rutile, which, if equilibrated with the magmas, are sufficient to impart the arc-type geochemical signature . Second, both northern and southern Lesser Caucasus magmas appear to be produced at significantly lower temperatures than the low K/Nb magmas of Eastern Anatolia, which are thought to be derived from the convecting mantle (McNab et al., 2018), such that the Lesser Caucasus magmas either have re-equilibrated in colder lithosphere or are entirely derived from the lithosphere.
Based on the smooth trend in 87 Sr/ 86 Sr versus 143 Nd/ 144 Nd in Fig. 11 and the gradational changes in trace element patterns (Fig. 8), it appears that magmas are produced from melting of two source types in varying proportions. If these two source types reflect geochemical end-members, one would be typified by the northern Lesser Caucasus volcanic rocks, whereas the other end-member would be most clearly seen in NW Iran.
Unfortunately, constraints on the lithospheric thickness are very limited for the northern Lesser Caucasus (Fig. 1), with the Priestley et al. (2012) model showing only that it is <100 km. Taking a typical conductive geothermal gradient for this lithosphere, melting would be expected to occur close to the lithosphere-asthenosphere boundary. Mantle-melt equilibration close to the base of the crust (Fig. 17) suggests that the lithosphere is very thin. If this were the case it is likely that there would be some melting of the convecting mantle based on elevated mantle T p in the region (McNab et al., 2018). However, such a thin lithosphere would make it difficult to insulate the crust from heating, such that more crustal contamination would be expected. As was noted above, this mantle source type appears also to be present in the southern Lesser Caucasus, where the lithosphere is sufficiently thick to suppress melting of anhydrous convecting mantle, but where melting of a hydrous lithospheric source would still be possible. It is likely that this magma type is derived from melting of the base of the lithosphere in response to small-scale convective removal, as suggested by Neill et al. (2015). This is shown in Fig. 18a as melt zone 1, with small portions of lithosphere being delaminated. This type of delamination is suggested to occur because a very small amount of water (a few hundred ppm) left over from previous subduction lowers the viscosity of the mantle sufficiently to allow more vigorous convection to render the lithosphere-asthenosphere boundary unstable (Kaislaniemi et al., 2014). Small-scale convective removal is preferred over catastrophic large-scale delamination because volcanism in the Lesser Caucasus is generally small-scale and sustained, with no evidence of crustal contamination, consistent with a continuous, less invasive process. Although this magma type is a mixing component in the southern Lesser Caucasus, it is not seen in NW Iran because the lithosphere is too  (Hu et al., 2017). Data filtering follows the approach outlined by Hu et al. (2017). Each volcanic highland has a range of Moho depth estimates of around 20 km. For both the Sr/Y and (La/Yb) N the median value for each volcanic highland is taken and then the two values are averaged to give an estimate of crustal thickness. Lithosphere thickness estimates are from Priestley et al. (2012). It is noted that they only estimate thicknesses where the lithosphere is >100 km, so the lithospheric thickness in the NW is schematic after Neill et al. (2015). Melting at the wet peridotite solidus lower depth limit is on the basis of upwelling convecting mantle at 1300-1400 C T p , which would lead to the wet solidus being crossed at $140 km depth (Green & Falloon, 2005). Also shown are the two melting zones (labelled 1 and 2) that are discussed in the text. Melting zone 1 is along the base of the lithosphere, and is suggested to be in response to small-scale convective removal of the lithosphere. Melting zone 2 is in the mid-lithosphere. (b) Sketch of the lithospheric structure of NW Iran and the location of melting in the mid-lithosphere. (c) Schematic illustration of the thermal relaxation of a kinked geothermal gradient leading to melting in the mid-lithosphere from initial underplating at time t 0 , through thermal evolution at time t 1 , to establishing a cratonic geotherm by time t 2 . thick for wet or dry melting at its base (Green & Falloon, 2005;Priestley et al., 2012). As the lithosphere is thickened, the degree of melting in response to convective removal will decrease as the wet peridotite solidus is approached with depth (Fig. 18a), which may help to explain the lower degrees of melting in the south.
The second melting zone is the sole magma source for volcanic rocks in NW Iran (Fig. 18b). Despite very different lithosphere thicknesses, both NW Iran and southern Lesser Caucasus magmas formed under similar conditions within the lithosphere (Fig. 17). It is suggested that continental collision could have led to the formation of a kinked geothermal gradient (Fig. 18c). A linear craton-style geothermal gradient is likely to be unrealistic for orogenic lithosphere . It has been suggested that orogenesis may proceed by underthrusting of mantle lithosphere from the oncoming plate (e.g. Willett et al., 1993). In this case underthrusting of Arabian mantle lithosphere, with its own pre-existing geothermal gradient, may lead to a kinked geothermal gradient with the kink forming along the subduction plane between overlying Eurasian lithosphere and underlying Arabian lithosphere. Upon initial collision, this kink will be very sharp (t 0 in Fig. 18b), but over time the gradient will thermally relax (t 1 ) towards a gradient that resembles the one characteristic for cratons (t 2 ). This thermal relaxation would heat the lithosphere just under the kink, which would previously have been cool. Dewatering and melting of this horizon would follow, which is what we argue we see in our P-T estimates of magma generation in Fig. 17. Given the lower degrees of melting in the southern Lesser Caucasus, it appears that this mechanism develops only low-degree melts.
To summarize, in the northern Lesser Caucasus magma is generated by small-scale delamination events heating the base of the lithosphere. This process also occurs in the southern Lesser Caucasus, but at greater depths. Magmas generated by this mechanism mix with a second type of magma produced in the midlithosphere from thermal relaxation of a kinked geotherm. Further south, in NW Iran, melting of the base of the lithosphere is suppressed owing to the depth being too great even for wet melting. However, melting in the mid-lithosphere continues.
How does thicker mantle lithosphere influence the composition of the mantle source?
The gradient in Sr-Nd isotope ratios in Fig. 11, along with a lack of evidence for crustal contamination, requires that there be some variation in the composition of the mantle source. One possibility is that magmas are tapping different lithosphere domains. Northern Lesser Caucasus magmas are exclusively found on Mesozoic arc lithosphere (Fig. 2), whereas southern Lesser Caucasus magmas are on South Armenian Block lithosphere, or else are very close to the suture. However, crossing such lithospheric sutures would probably result in a step-change in isotope ratios, rather than the smooth gradation observed in the Lesser Caucasus. A minority of samples from the Vardenis volcanic highland (Figs 1 and 2) are thought to be on the east (Mesozoic arc) side of the suture; these are shown by an asterisk in Table 3. As Table 3 shows, these samples have Sr-Nd isotope compositions indistinguishable from those of other Vardenis samples, suggesting that the suture zone is not a major dividing line in isotope composition. This suggests that melting zone 2 in the mid-lithosphere has a different composition from melting zone 1 at the base of the lithosphere. In this discussion on the nature of the mantle sources, only the most mafic samples (>4 wt % MgO, <54% SiO 2 ) are used to try and minimize the effects of fractional crystallization on trace element contents.
The lower crust is thought to behave as a weak layer during continental collision (Bü rgmann & Dresen, 2008). This could lead to some lower crust being incorporated into the mantle lithosphere perhaps during underthrusting of Arabian lithosphere, as has been suggested for numerical models of other collision zones such as the Himalayas (e.g. Toussaint et al., 2004;Li et al., 2011). This could lead to the addition of lower continental crust (LCC) to melting zone 2 in the middle of the lithosphere, and hence enrichment of the mantle source of the southern Lesser Caucasus. Incorporation of LCC into the mantle source could result in significant enrichments in all of the most incompatible trace elements, as in Fig. 8e, given the much higher concentrations of most incompatible trace elements in the lower crust relative to the primitive mantle (Sun & McDonough, 1989;Rudnick & Gao, 2003).
The La/Nb ratios of the least enriched samples from the northern Lesser Caucasus and the same ratio estimated for the LCC (Rudnick & Gao, 2003) provide a Lower continental crust (LCC) has a La/Nb ratio with a similar value of around 1Á6 (Rudnick & Gao, 2003). Therefore if LCC is enriching the melt source, the La/Nb ratio should be close to invariant. serendipitous coincidence, with both around 1Á6 (Fig. 19). Given that La/Nb is unaffected by melting and crystallization (i.e. both elements have similar bulk partition coefficients with respect to the melt) regardless of source mineralogy (Ionov et al., 2002;McKenzie & O'Nions, 1991), any variations in the La/Nb ratio should reflect variations in the source La/Nb. If the lower crust was responsible for the source enrichment, La/Nb should be near-constant along the transect. Samples from the southern Lesser Caucasus have much higher La/Nb ranging from 2Á5 to 5 (Fig. 19), which cannot be explained by the addition of average LCC.
However, it is worth noting that several estimates of the composition of individual lower crust sections do give higher La/Nb of up to 4Á9 (Weaver & Tarney, 1984;Villaseca et al., 1999;Jagoutz & Schmidt, 2013), such that the involvement of lower crust in the southern Lesser Caucasus magma source cannot be precluded.
Lithospheric mantle that is significantly shallower than the lithosphere-asthenosphere boundary (melting zone 2) is expected to be colder than the deep lithosphere (melting zone 1) prior to any heating, even if the lithosphere does have a kinked geothermal gradient (Fig. 17). These lower temperatures should stabilize minerals such as amphibole and phlogopite, which can retain chemical components derived from mantle metasomatism (Luth, 2003;Frost, 2006;and references therein). In the deep lithosphere, these components would be subject to upward mobilization by fluid release following dehydration at higher pre-melting temperatures. Underthrusting of Arabian lithosphere, as suggested in Fig. 18, would also introduce a new lithospheric domain, which could have a different composition (including isotopically) from the Eurasian lithosphere.
If this metasomatic material is responsible for the enrichments in the southern Lesser Caucasus, then it must have a composition capable of producing those most enriched melts. La/Yb and Sr/Y are two ratios that increase most dramatically with latitude, as shown in Fig. 20a and b. Both of these ratios show excellent correlations with P 2 O 5 content ( Fig. 20c and d) within all Lesser Caucasus samples. The La/Yb ratios and P 2 O 5 contents of the NW Iran samples (interpreted to be pure melting zone 2 magmas) appear as end-members on the Lesser Caucasus mixing line (Fig. 20d), suggesting that apatite may be a metasomatic phase in the host rock of melting zone 2 (Fig. 18). This is consistent with the presence of large apatite crystals in some southern Lesser Caucasus samples (up to 2 mm across; Fig. 5b). It is also consistent with previous suggestions that metasomatic apatite may be an important phase in the Iranian subcontinental lithospheric mantle (Pang et al., 2013). As was shown in Fig. 15, it is difficult to produce such high LREE concentrations (shown by La in Fig. 15) simply by lower degrees of melting or more garnet in the mantle source. Apatite is soluble in melts with a low SiO 2 content (Watson, 1980), such that apatite may simply be added to the initial modelled melt composition, as shown by the blue dashed line in Fig. 15, producing a much improved agreement with natural southern Lesser Caucasus samples.
Although apatite can explain the enrichment in LREE, it cannot explain the large enrichment in other elements such as Ba, which probably requires wider metasomatic reworking. Two minerals that could explain this enrichment in Ba are amphibole and phlogopite. As Fig. 21 shows, the high Ba/Rb ratio favours amphibole as the major metasomatic phase in the source of both the Lesser Caucasus and NW Iran magmas. However, in both the southern Lesser Caucasus and NW Iran there are a minority of samples for which phlogopite could be an important metasomatic phase. It is noteworthy that the presence of phlogopite in the mantle source is indicated in a minority of cases for Vardenis, but not Syunik (Fig. 21), given that biotite is found as a phenocryst phase only in mafic samples from Vardenis (see Petrography section). The presence of these metasomatic signatures in the southern Lesser Caucasus and NW Iran is consistent with melting zone 2 being host to apatite, amphibole, and occasionally phlogopite prior to melting.
Sr/Y is a ratio that does not follow a simple mixing pattern. It increases across the Lesser Caucasus, but is actually lower in NW Iran compared with the southern Lesser Caucasus. This ratio probably reflects the involvement of garnet in the mantle source (Defant & Drummond, 1990), which should become dominant at greater depths. In NW Iran, all melting is occurring in the mid-lithosphere; however, in the southern Lesser Caucasus magmatism is driven by melting in both the mid-lithosphere and at the base of the lithosphere (Fig. 18). This means that in the southern Lesser Caucasus, the average depth of melting may actually be greater than in NW Iran, despite the thinner lithosphere, leading to the mantle source being more dominated by garnet, increasing the Sr/Y ratio of the resulting magmas.
High Sr/Y ratios have been associated with adakites and the melting of oceanic slabs (Defant & Drummond, 1990). Adakitic signatures have been seen in some volcanic suites from NW Iran (e.g. Ghalamghash et al., 2016), suggesting that the Iranian enriched signature could somehow derive from slab melting. These adakites are generally andesites and dacites. The more primitive NW Iran lavas of Allen et al. (2013) and the southern Lesser Caucasus lavas have higher Y contents than true adakites, and it is likely that the adakite-like compositions of andesites and dacites are derived from fractional crystallization processes (e.g. Chiaradia et al., 2009), rather than slab melting.
Melting zone 2 does appear to have a more enriched composition than melting zone 1, as demonstrated by signatures of the metasomatic minerals apatite and amphibole in magmas derived from the southern Lesser Caucasus and NW Iran. However, variations between the chemistry of the northern and southern Lesser Caucasus volcanic rocks do not just reflect mixing of source reservoirs, but also changes in the degree of melting (lower degree of melting will increase La/Yb and Sr/Y ratios), and the mineralogy of the melt source (garnet in the melt source of the southern Lesser Caucasus also increases Sr/Y and La/Yb). The multiple parameters controlling the trace element composition of Lesser Caucasus volcanic rocks mean that variations between the northern Lesser Caucasus, southern Lesser Caucasus and NW Iran are often non-linear.

SUMMARY AND CONCLUSIONS
Magmas generated in the thicker lithosphere of the southern Lesser Caucasus have higher incompatible trace element concentrations, higher 87 Sr/ 86 Sr ratios and lower 143 Nd/ 144 Nd ratios than volcanic rocks from the northern Lesser Caucasus. A lack of consistent variation between the isotope compositions of basalts and rhyolites and SiO 2 suggests that crustal contamination is unimportant in generating the enriched geochemistry. The negative Nb-Ta anomalies, and enrichments in large ion lithophile elements and LREE are instead likely to be produced by partial melting of a subductionmodified mantle source. This subduction component is uniform across the Lesser Caucasus and is probably inherited from Mesozoic (Tethyan) slab subduction prior to continental collision. The more enriched geochemistry of southern Lesser Caucasus rocks is the result of lower degrees of melting, an increased proportion of garnet in the mantle source, and also a distinct mantle source composition.
The temperatures of melt formation in the mantle are all less than 1200 C, which when compared with the much higher temperatures for magmas formed in the asthenosphere below nearby Eastern Anatolia suggests that magma generation occurs in the lithosphere, which is also consistent with the ubiquitous subduction signature. Very similar conditions of melt generation in NW Iran and the southern Lesser Caucasus, as well as several similarities in geochemistry, suggest a common magma generation mechanism in the mid-lithosphere, despite very different lithospheric thicknesses. This magma type appears to mix in the southern Lesser Caucasus with magma from a second source, which probably originates at the base of the lithosphere. This latter magma source is the sole site of magma generation in the northern Lesser Caucasus, where melting occurs at the base of an $50 km thick lithosphere. In the southern Lesser Caucasus melting occurs at $75 km depth, significantly shallower than the estimated 120 km thickness of the lithosphere. Melting in the midlithosphere occurs as a result of relaxation of a kinked geothermal gradient, whereas melting at the base of the lithosphere is the result of small-scale delamination events. This latter melting mechanism proceeds only until the point where the lithosphere becomes too thick to melt at its base, even if the mantle peridotite is hydrated. The enriched composition of the midlithospheric mantle source could be derived from the incorporation of weak lower crust during collision. However, several signatures of the metasomatic minerals amphibole, apatite and occasionally phlogopite suggest that the enriched nature of the mantle source in the mid-lithosphere is derived from the increased retention of metasomatic components in hydrous minerals prior to the post-collisional magmatism.
Interestingly, it appears that a melt source exists at somewhat less than 100 km depth regardless of the lithospheric thickness across the Anatolian-Armenian-Iranian plateau. This is consistent with the geophysical observations of Maggi & Priestley (2005), which show a low shear wave velocity at around 100 km depth below the entire plateau. Further work on understanding the interplay between lithospheric thickness and melt generation in continental collision zones would benefit from detailed tomographic work in the critical region of the Lesser Caucasus to help us better understand how the thickness of the lithosphere varies along this mountain range. Investigations of the petrogenesis of primitive magmatic rocks from NW Iran could elucidate whether thermal relaxation of a kinked geothermal gradient is a viable mechanism to generate magma in the mid-lithosphere. Studies of stable fluid-sensitive isotope systems such as O and B would help decipher the nature and role of inherited subduction components in the generation of collision-related magmas.