The Canary Islands hot spot: New insights from 3D coupled geophysical–petrological modelling of the lithosphere and uppermost mantle

Earth and Planetary Science Letters 409 (2015) 71–88

Contents lists available at ScienceDirect

Earth and Planetary Science Letters www.elsevier.com/locate/epsl

The Canary Islands hot spot: New insights from 3D coupled geophysical–petrological modelling of the lithosphere and uppermost mantle Javier Fullea a,∗ , Antonio G. Camacho a , Ana M. Negredo a,b , José Fernández a a b

Institute of Geosciences (CSIC, UCM), Plaza de Ciencias, 3, ES-28040 Madrid, Spain Dept. of Geophysics, Facultad CC. Físicas, Universidad Complutense de Madrid, Plaza de Ciencias, 1, ES-28040 Madrid, Spain

a r t i c l e

i n f o

Article history: Received 23 June 2014 Received in revised form 9 October 2014 Accepted 15 October 2014 Available online xxxx Editor: J. Brodholt Keywords: geophysical–petrological modelling mantle composition seismic tomography Canary Islands mantle plume hot spot

a b s t r a c t The Canary archipelago (NW Atlantic African margin) is one of the best studied volcanic chains in the world yet its structure and geodynamic evolution are still under considerable debate. Oceanic island volcanoes typically form over hot spots due to upwelling of plume material followed by decompression melting and melt migration up to the surface. Here, the 3D lithospheric-uppermost mantle thermochemical structure beneath the Canary Islands is studied using an integrated and selfconsistent geophysical–petrological approach exploiting the wealth of available data after decades of geophysical and petrological studies plus recent satellite data. A precise knowledge of the presentday thermal and compositional mantle structure beneath the Canary Islands is a key element to understand the geodynamic evolution of the area and, on a global scale, the thermal state of the Earth’s mantle beneath hot spots. Our results suggest a likely chemically depleted and mechanically strong lithosphere showing no significant thinning with respect to the surrounding oceanic and continental domains (110 ± 20 km thick). Models without a positive temperature anomaly in the sub-lithosphere (characterized by mantle T pot = 1335 ◦ C) fail to reproduce the observed sub-lithospheric seismic anomaly over the Canary Islands. A thermal sub-lithospheric anomaly of +100 ◦ C (mantle potential temperature of 1435 ◦ C) with respect to ambient mantle beneath the Canaries is able to explain both observed seismic tomography anomalies and measured geophysical and geodetic data. Such a sub-lithospheric thermal anomaly requires a dynamic contribution of 150–400 m to the static topography to match the presentday observed elevation in the Canary Islands and associated swell. © 2014 Elsevier B.V. All rights reserved.

1. Introduction Oceanic hot spots are generally related to deep thermal anomalies of +(100–300) ◦ C with respect to ambient mantle temperatures regardless of possible additional bulk chemical changes or volatile enrichment (e.g., Herzberg and Asimow, 2008; White, 2010). Oceanic island volcanoes typically form over hot spots due to upwelling of plume material followed by decompression melting and melt migration up to the surface. The Canary Islands (NW Atlantic African margin) are one of the best studied volcanic chains in the world yet some aspects of its structure and geodynamic evolution are still controversial (e.g., Anguita and Hernán, 2000). The Canary archipelago shows some specific characteristics with respect to other mid oceanic volcanic chains classically explained

*

Corresponding author. E-mail address: [email protected] (J. Fullea).

http://dx.doi.org/10.1016/j.epsl.2014.10.038 0012-821X/© 2014 Elsevier B.V. All rights reserved.

by the mantle plume hypothesis (e.g., Hawaii): (i) lack of a prominent bathymetric swell (e.g., Watts, 1994); (ii) long term and irregular volcanic evolution of >20 Ma and, perhaps 70–80 Ma in Fuerteventura (Le Bas et al., 1986); (iii) low melt production rates (Hoernle and Schmincke, 1993) and multiple cycles of volcanic activity. Alternative hypotheses on the origin of the Canary Islands excluding a mantle plume (see Anguita and Hernán, 2000 for an overview and references) are a propagating fracture connecting the archipelago and African Atlas Mountains, compression-related tectonic uplift, and local rifting in the Islands. Interaction between a plume and small-scale edge driven convection has also been suggested (e.g., Geldmacher et al., 2005). Most of the proposed geodynamic scenarios integrate as a key element the existence of a deep thermal sub-lithospheric anomaly (i.e., mantle plume/broad thermal anomaly) and its interaction with old and slowly moving Jurassic oceanic lithosphere close to the north Atlantic African passive margin. In this work we analyze the lithospheric-uppermost mantle thermochemical structure beneath the Canary Islands using

72

J. Fullea et al. / Earth and Planetary Science Letters 409 (2015) 71–88

eruptive process took place offshore, October 2011–March 2012, near the southern shoreline of El Hierro (González et al., 2013). 3. Geophysical and petrological setting Geophysical studies in the Canary archipelago at a crustal scale include wide-angle deep seismic experiments (e.g. Banda et al., 1981; Ye et al., 1999; Dañobeitia and Canales, 2000; see locations in Fig. 1), seismic tomography (Krastel and Schmincke, 2002; García-Yeguas et al., 2012) and receiver functions (Lodge et al., 2012; Martinez-Arevalo et al., 2013), magnetotellurics (Pous et al., 2002), gravity modeling (e.g., Ranero et al., 1995; Montesinos et al., 2006; Camacho et al., 2009, 2011), and elastic thickness estimates (e.g. Watts, 1994; Watts et al., 1997; Dañobeitia et al., 1994; Canales and Dañobeitia, 1998). Fig. 1. Elevation map (Smith and Sandwell, 1994; Smith and Sandwell, 1997) of the study area. Red lines show the location of published seismic lines used to constrain the crustal structure (see Appendix A). a1 and a2: CD82-P11 and CD82-P12 respectively (Watts et al., 1997; Dañobeitia and Canales, 2000). a3: M24-P1 (Ye et al., 1999). a4: Banda et al. (1981); Dañobeitia and Canales (2000). Lithospheric mantle compositional domains (black dotted lines) based on crustal tectonics, and petrological and geophysical considerations are shown (see text for further details). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

an integrated geophysical–petrological approach, “LitMod” (Afonso et al., 2008; Fullea et al., 2009), able to reduce the uncertainties associated with the modelling of different data sets separately, avoid inconsistencies, and exploit the different sensitivities of geophysical observables. The wealth of available data after decades of geophysical and petrological studies (gravity and geoid anomalies, elevation, seismic and mantle xenoliths) is exploited here, along with recently released satellite data (GOCE gravity gradients), to constrain the present-day thermal and compositional 3D structure of the lithosphere/uppermost mantle beneath the Canary Islands as a fundamental element to understand its geodynamic evolution and, on a global scale, the thermal state of the Earth’s mantle beneath hot spots.

3.1. Crustal seismic structure The Canary Islands were formed over Jurassic oceanic crust, characterized by a thickness of 5–7 km (Banda et al., 1981), which progressively thickens towards the Atlantic Moroccan passive margin to 27–35 km (Contrucci et al., 2004; Klingelhoefer et al., 2009; Spieker et al., 2014). The seismic structure beneath the central islands is defined by a volcanic edifice (V p = 5.5–6 km/s) of variable thickness underlain by a 6–7-km-thick lower crust (V p = 6.6–7.3 km/s) (Ye et al., 1999; Dañobeitia and Canales, 2000). In the eastern islands the upper crust (V p = 6–6.7 km/s) is 5–8 km thick and the lower crust is absent (Dañobeitia and Canales, 2000 and references therein). Furthermore, under the central and eastern islands a 7–12-km-thick layer defined by P-wave velocities (7.4–8 km/s), higher than those of typical lower crust but lower that the average Jurassic oceanic uppermost mantle velocities in the neighborhood of the islands (8 km/s), has been interpreted as a magmatic underplating (Dañobeitia and Canales, 2000; Freundt and Schmincke, 1995; Lodge et al., 2012) (see Appendix A for more details). 3.2. Lithosphere and uppermost mantle structure and composition

2. Geological setting The Canary archipelago, limited to the west and east by magnetic anomalies M25 and S1 (Verhoef et al., 1991; Roest et al., 1992), lies on Jurassic (150–170 Ma) oceanic lithosphere adjacent to the NW African passive margin (Fig. 1). This margin hosts a 3000-km-long volcanic belt that includes a considerable number of seamounts and volcanic islands. The seven major islands in the Canary archipelago exhibit a long volcanic history (70–80 Ma) and hence multiple oceanic volcanic islands stages (i.e., seamount, shield, erosional) are well represented. The volcanic activity in the Canary archipelago shows East–West age progression with the oldest exposed volcanic rocks in Fuerteventura (20 Ma, up to 70–80 Ma according to Le Bas et al., 1986) and the youngest (<4 Ma) in the western islands (La Palma and El Hierro). The eastern islands, Fuerteventura and Lanzarote, are parallel to the NW African margin and show a rather flat topography (max elevation of 807 and 607 m respectively) characteristic of its erosional stage. These two islands along with the conception Bank, north of Lanzarote, define the East Canary Ridge (Ancochea et al., 2004 and references therein). The central islands, Gran Canaria, La Gomera and Tenerife, exhibit an E–W trend. Tenerife and Gran Canaria are in the post-shield stage with rejuvenated volcanism which is absent in La Gomera (erosional stage, no volcanism in the last 2–3 Ma). The young western islands, La Palma and El Hierro, align along an N–S trend and are currently at a rather juvenile shield stage (Ancochea et al., 2004 and references therein). The most recent

Lithospheric-uppermost mantle scale studies in the Canary Islands are relatively scarce. Based on 2D seismic reflection and gravity data, Ranero et al. (1995) modelled a moderate lithospheric thinning from a lithosphere–asthenosphere-boundary (LAB) depth of 100 km west of the Canaries in the Jurassic oceanic lithosphere to about 80 km in the western islands (La Palma and El Hierro). These authors argued that the topographic swell associated with the archipelago could not be explained by crustal variations and that deep density anomalies were required instead. Neumann et al. (1995) estimated a lithosphere thickness of only 27 km under Lanzarote based on petrological constraints, and suggested thermal erosion as a possible explanation. A recent seismic tomography model based on multimode inversion of surface- and S-wave forms in Europe shows a low velocity anomaly area centered in the central islands and affecting the whole Canarian domain at lithospheric (50–110 km) and sublithospheric depths (150–260 km) (Legendre et al., 2012). Earlier body-wave seismic tomography models have also identified a deep and broad low velocity zone in the mantle beneath the Canary Islands (Hoernle et al., 1995). A recent multiple-frequency P-wave velocity tomography model also shows a negative velocity anomaly underlying the Canary Islands in the lithosphere and upper mantle (Bonnin et al., 2014). The negative anomaly in Bonnin et al. (2014) seismic model seems to be restricted to the western islands in the lithosphere, progressively shifting north-westwards of La Palma in the mantle transition zone (their Fig. 6). However, the extent to which this W–E mantle velocity pattern in the upper mantle is a

J. Fullea et al. / Earth and Planetary Science Letters 409 (2015) 71–88

73

Table 1 Bulk mantle compositions used in this work.

SiO2 Al2 O3 FeO MgO CaO Na2 O Mg#

(1) OI1-u: Av. Harz. Lanzarote (wt%)a

(2) OI1-u: Av. HEXO Tenerife (wt%)b

(3) OI1-u: Av. Harz. La Palma av. (wt%)c

(4) OI2-u: Av. HLCO Tenerife (wt%)b

(5) Av. Ocean floor peridot. (wt%)d

(6) Av. Azrou Middle Atlas (wt%)e

(6) PUM M&S95 (wt%)f

43.78 0.7 7.79 46.1 0.6 0.1 91.34

43.32 0.61 8.04 45.31 0.81 0.14 90.96

43.07 0.53 8.43 45.19 0.68 0.17 90.53

42.14 0.73 8 .8 44.14 1.68 0.18 89.94

45.09 2.33 8 .4 41.23 1.32 0.23 89.7

43.48 2.38 8 42.6 2.83 0.24 90.47

45 4.45 8.05 37.8 3.55 0.36 89.3

a

Average composition of Lanzarote spinel harzburgite suite from Neumann et al. (1995). Average composition of xenoliths from Tenerife. HEXO harzburgites show strongly exsolved and deformed opx porphyroclasts, and HLCO harzburgites and lherzolites contain clear opx porphyroclasts (Neumann et al., 2002). c Xenolith data averaged from La Palma in Wulff-Pedersen et al. (1996). d Average major element (Shibata and Thompson, 1986) and modal (Dick, 1989) compositions of ocean floor peridotites. e Average composition from Azrou Middle Atlas xenoliths from Wittig et al. (2010). f PUM stands for Primitive Upper Mantle, M&S95 refers to McDonough and Sun (1995). In (1–4) OI1-u and OI2-u correspond to ultra-refractory mantle compositions in oceanic islands with and without petrographic signs of metasomatism (either in situ in the mantle or resulting from interaction with the host magma during ascent) respectively (Simon et al., 2008). b

well resolved feature is still controversial due to the high level of noise in the seismic stations over the Canary Islands (A. Villaseñor, personal communication). The mantle composition sampled from mantle xenoliths erupted in the Canaries reflects a highly refractory mantle affected by significant melt extraction and, in some cases, metasomatism (e.g., Neumann et al., 2002; Simon et al., 2008). There are no mantle xenoliths directly sampling the lithosphere from the neighbor Moroccan margin in our study region. Farther to the north east, in the Moroccan Middle Atlas, there are abundant mantle xenoliths associated with Cenozoic volcanic fields (Wittig et al., 2010; Raffone et al., 2009). These xenoliths are coarse-grained spinel perdiotites with no signs of serpentinization, generally lherzolitic with some exceptional harzburgites. In contrast to the mantle beneath the islands, the lithosphere from the Middle Atlas xenoliths appears to be relatively fertile (Wittig et al., 2010; Raffone et al., 2009) and close to global average composition for Phanerozoic lithospheres (e.g., Griffin et al., 1999) (Table 1). Xenolith samples from the Canary Islands are generally spinel harzburgites and rarely spinel lherzolites and dunites (e.g. Wulff-Pedersen et al., 1996). The modal olivine–orthopyroxene–clinopyroxene ratios, mineral and whole rock chemistry are similar in the harzburgite and lherzolite xenoliths in La Palma, El Hierro and Lanzarote, reflecting a highly depleted mantle (Neumann et al., 1991; Neumann et al., 1995; Wulff-Pedersen et al., 1996). Xenoliths from Tenerife appear to have undergone more intense metasomatism compared to the other islands in the archipelago (e.g., Neumann et al., 1995, 2002; Wulff-Pedersen et al., 1996). Spinel lherzolites and harzburgites from Tenerife present two types of textures related to the appearance of opx: (i) HEXO harzburgites show strongly exsolved and deformed opx porphyroclasts, and (ii) HLCO harzburgites and lherzolites containing clear opx porphyroclasts (Neumann et al., 2002). HLCO xenoliths contain more ol and cpx and less opx than HEXO xenoliths which are similar in bulk chemistry to xenoliths from La Palma, El Hierro and Lanzarote: HLCO xenoliths are relatively CaO rich with respect to HEXO and the xenoliths from the other Canary Islands (1.7 wt% vs 0.6–0.8 wt% respectively) (Table 1). Major elements differences between Tenerife HEXO (and La Palma, El Hierro and Lanzarote) and HLCO xenoliths can be seen in the light of the general classification for ocean islands given by Simon et al. (2008): Tenerife HEXO falls within OI1-u group whereas Tenerife HLCO corresponds to OI2-u group. OI2-u and OI1-u are ultra-refractory mantle compositions in oceanic islands with and without petrographic signs of metasomatism (either in situ or resulting from interaction with the host magma during ascent) respectively. Hence, OI1-u represents the pre-ocean-island

composition before the onset of the ocean island magmatism (Simon et al., 2008). With respect to a typical ocean floor peridotitic composition, the xenoliths from the Canaries are all depleted in Al2 O3 (2.33 wt% vs 0.5–0.7 wt% respectively). In contrast, xenoliths from the continental Middle Atlas are CaO and Al2 O3 rich with respect to Canarian xenoliths (2.83 wt% and 2.38 wt% vs 0.6–0.8 wt% and 0.5–0.7wt % respectively) (Table 1). 3.3. Effective elastic thickness and lithospheric flexure Geophysical studies have led to diverging conclusions regarding the mechanical strength of the lithosphere beneath the Canary Islands. Watts (1994) and Watts et al. (1997) found, on the basis of seismic and free air anomaly data, a relatively weak (elastic thickness T e = 20 km) lithosphere beneath the archipelago compared to the expected value for unperturbed Jurassic oceanic lithosphere (e.g., T e = 50 km, Filmer and McNutt, 1989), probably related to thermal weakening produced by thermal erosion. A similarly low value (T e = 23 km) was determined by Dañobeitia et al. (1994) assuming an isotropic elastic plate. However, these authors found that the neighboring African passive margin influences the mechanical behavior of the Canarian lithosphere (i.e., regional NW–SE flexural trend), leading to a higher estimate of T e = 35 km. Based on the spectral relationship between topography and gravity Canales and Dañobeitia (1998) calculated a T e value of 28–36 km. T e estimates given by Canales and Dañobeitia (1998), Filmer and McNutt (1989) and Dañobeitia et al. (1994) (T e = 30–50 km) would fall within the range of expected values for unperturbed (150–170)-Ma-old oceanic lithosphere (e.g. Watts, 2001). In contrast, a value as low as T e = 20 km (Watts, 1994; Watts et al., 1997) would imply thermal rejuvenation of the lithospheric mantle. 4. Geophysical and geodetic observables: constraining data Regional geophysical and geodetic observables used as constraining data in this work were taken from different sources (Fig. 2). Geoid and free-air anomalies were obtained from the global Earth model EGM2008 (Fig. 2A and B) (Pavlis et al., 2008). The geoid signal has been filtered to remove long wavelengths of deep origin (>4000 km, degrees 2–9) and retain the effects of density anomalies shallower than ∼400 km depth (Bowin, 2000). Bouguer anomalies in the islands came from different sources: La Palma (Camacho et al., 2009), Tenerife (Vieira et al., 1986; Ablay and Kearey, 2000; Araña et al., 2000; Gottsmann et al., 2008), Gran Canaria (Anguita et al., 1991; Camacho et al., 2000),

74

J. Fullea et al. / Earth and Planetary Science Letters 409 (2015) 71–88

Fig. 2. Geophysical and geodetic observables. (A) Geoid anomaly from the global Earth model EGM2008 (Pavlis et al., 2008). Long wavelengths (n < 9) have been removed to retain the lithospheric signal. (B) Free air anomaly from satellite data (EGM2008, Pavlis et al., 2008). (C) Bouguer anomalies in the islands and Morocco come from different sources (see text). For the rest of the African continental areas and offshore, Bouguer anomaly is computed applying the complete topographic correction to the to free-air satellite data (EGM2008) assuming a reduction density 2200 kg/m3 . (D) Elevation from ETOPO2 DEM (Smith and Sandwell, 1994; Smith and Sandwell, 1997). (E1–E6) Gravity gradients from the satellite-only (mission GOCE) global Earth model GOCO03S (http://www.goco.eu/), computed up to degree and order 220 at 255 km above the sea level (GOCE mission perigee height) using a spherical harmonics synthesis code (see appendix in Fullea et al., 2015 and Martinec and Fullea, 2015). The sub-indexes of V denote second derivative of the Earth’s potential (i.e., gravity gradient) along the south–north, east–west, and vertical directions (x, y, z respectively).

J. Fullea et al. / Earth and Planetary Science Letters 409 (2015) 71–88

Lanzarote (Camacho et al., 2001) and Fuerteventura (Montesinos, 1999). In Morocco, Bouguer anomalies come from Hildenbrand et al. (1988). For the rest of the African continental areas and offshore, Bouguer anomaly were computed from free-air satellite data (EGM2008) corrected using the software FA2BOUG for a reduction density of 2200 kg/m3 (Fullea et al., 2008) (Fig. 2C). The elevation data, i.e., topography and bathymetry, came from the ETOPO2 Global Data Base (Fig. 2D) (Smith and Sandwell, 1994, 1997). Gravity gradients were taken from the recent global Earth model GOCO03S (http://www.goco.eu/) based on GOCE satelliteonly data, computed up to degree and order 220 (lateral resolution of about 90 km) at the satellite height (≈255 km) using a spherical harmonics synthesis code (see appendix in Fullea et al., 2015 and Martinec and Fullea, 2015) (Fig. 2E1–E6). 5. The method: integrated geophysical–petrological modeling A general description of the method used in the present study (LitMod approach) is given by Afonso et al. (2008) and Fullea et al. (2009). Here we present a general overview of the fundamentals for completeness, with a special focus on the topics relevant to our study. 5.1. The lithosphere–asthenosphere boundary In this paper we adopt the definition of the lithosphere– asthenosphere boundary (LAB) based primarily on the temperature and compositional distributions. Therefore, we assume that the lithospheric mantle is defined: (i) thermally, as the portion of the mantle characterized by a conductive geotherm, and (ii) compositionally, as the portion of the mantle characterized by a, generally, different (normally, more depleted) composition with respect to the fertile primary composition in the sub-lithosphere (i.e., PUM in Table 1). 5.2. The geotherm The lithospheric geotherm is computed under the assumption of steady-state heat transfer in the lithospheric mantle, considering a P–T-dependent thermal conductivity (Afonso et al., 2008; Fullea et al., 2009). Between the lithosphere, with its base defined by the 1315 ◦ C isotherm, and sub-lithosphere mantle a “transition” region (buffer layer) with variable thickness and a continuous linear super adiabatic gradient is assumed (i.e., heat transfer is controlled by both conduction and convection processes, see Fullea et al., 2009 for details). Below the buffer layer the geotherm is described by an adiabatic temperature gradient forced to be in the range 0.35–0.6 ◦ C/km. The mantle potential temperature, T pot , is defined as the temperature that the convecting mantle would attain if it were adiabatically decompressed to the surface without melting (McKenzie and Bickle, 1988). The parameterization of the temperature distribution at the LAB and buffer layer assumed here (see Fullea et al., 2009 for details) corresponds to a T pot ≈ 1335 ◦ C in the adiabatic mantle.

75

5.4. Forward modelling: geophysical and geodetic observables A detailed description on the calculation of synthetic gravity and geoid anomalies, and gravity gradients for a given 3D density distribution can be found in Fullea et al. (2009, 2015). The predicted/synthetic surface elevation in each model column (its buoyancy) is determined according to local isostasy by integrating the crustal and mantle densities from the surface down to the base of the model (400 km depth) and comparing subsequently with a calibration column (details on the calibration procedure are given in Afonso et al., 2008 and Fullea et al., 2009). A regional or flexural isostatic balance is computed using the finite-difference code TISC (García-Castellanos, 2002), which requires a 2D distribution of loads and a representative effective elastic thickness (T e ) for the lithosphere. In our approach the flexural loads are determined by pressure variations of lithospheric origin at the base of the model domain (see Appendix C2 in Fullea et al., 2009). Therefore, in the ideal case where the observed topography and the modelled lithospheric structure are in complete local isostatic equilibrium, the pressure at the compensation level would be constant throughout the entire model and the load distribution would be zero. Departures from local isostasy, i.e. lateral variations of the pressure at the bottom of the model or mismatches between observed and synthetic local isostasy elevation, are absorbed to a varying extent by the strength of the lithosphere depending on their wavelength and assumed T e value. 6. Modelling results The present study is focused on the regional scale thermochemical structure of the uppermost mantle beneath the Canary archipelago. Therefore, a reliable crustal structure based on previous studies needs to be defined and kept fixed throughout the modelling process (see Appendix A). Fixing the crustal structure allows us to explore a range of different thermal and compositional lithospheric mantle models, and to compare their outputs against measured geophysical and geodetic data (Section 4, Fig. 2) and other geophysical and petrological studies (i.e., seismic tomography, elastic thickness, mantle xenoliths, Sections 3.2 and 3.3). The crustal model considered in this study consists of five layers (characterized by different physical properties) (see Table 4): (1) sediments, (2) upper crust/volcanic edifice, (3) middle/oceanic crust, (4) lower crust and (5) magmatic underplating (Figs. 9 and 10). Further details on the modelled crustal structure are given in Appendix A. The lithospheric models presented here are characterized by two essential parameters: bulk chemical composition and temperature distribution (i.e., the geotherm). The geotherm, for a given crustal structure, is essentially determined by the depth of the LAB isotherm (see Section 5.2). In the following, we explore compositional and lithospheric thickness (i.e. LAB depth) variations below the Canary Islands, illustrated by a number of lithospheric models. 6.1. Mantle compositional domains

5.3. Thermodynamic framework Stable mineral assemblages in the mantle are calculated using a Gibbs free energy minimization as described by Connolly (2005). The composition is defined within the major oxide system NCFMAS (Na2 O–CaO–FeO–MgO–Al2 O3 –SiO2 ). All the stable assemblages in this study are computed using a modified/augmented version of the Holland and Powell’s (1998) thermodynamic database (Afonso and Zlotnik, 2011). The density and seismic velocities in the mantle are determined according to the elastic moduli and density of each end-member mineral as described by Connolly and Kerrick (2002) and Afonso et al. (2008).

Based on crustal tectonics, petrological and geophysical considerations, three broad lithospheric domains can be defined in the study region: (1) Oceanic lithosphere; (2) Canary Island lithosphere; and (3) Moroccan continental lithosphere (Fig. 1). Domain 1 comprises likely undisturbed Jurassic oceanic lithosphere; its bulk composition is taken from an average ocean floor peridotite (Table 1). The composition of domain 3 in this work is based on the average major elements chemistry from mantle xenoliths erupted in the Moroccan Middle Atlas (northeastwards of the study region) (Table 1). The xenoliths from Middle Atlas (Wittig et al., 2010; Raffone et al., 2009) are the closest Phanerozoic suite to the

76

J. Fullea et al. / Earth and Planetary Science Letters 409 (2015) 71–88

Fig. 3. Lithospheric vertical profiles of the seismic velocities and density associated with the mantle compositions discussed in this work. The temperature distribution corresponds with a 30-km-thick and 100-km-thick crust and lithosphere respectively. The solid black, red, blue, green, and dotted and dashed black lines correspond with lithospheric mantle compositions 1–6 in Table 1 respectively. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Atlantic Moroccan margin in our study region. Although some of the Middle Atlas samples show evidence of metasomatic alteration (Wittig et al., 2010), this is mostly reflected in minor and trace rather than major elements, the later controlling bulk density and seismic velocities relevant to this study. The major element composition from Middle Atlas used here (Table 1) is an average of eleven spinel-bearing perdotite samples (ten lherzolites and one harzburgite) from Azrou volcanic field (Wittig et al., 2010). This average composition from Azrou Middle Atlas suite is similar to average peridotitic compositions from global data for samples of similar tectonothermal age (e.g., Griffin et al., 1999). In the absence of direct mantle sampling of the Atlantic Moroccan margin adjacent to the Canary Islands, we take the average composition from Azrou Middle Atlas suite (Table 1) as a reasonable proxy to the mantle composition beneath our lithospheric domain 3 (Moroccan continental lithosphere) (Fig. 1). The Canary archipelago lies on domain 2, where the oceanic lithosphere seems to have been affected by significant processes of melt extraction and, to some extent, subsequent metasomatism (see Section 3.2 and Table 1). Mantle xenoliths available from the Canary Islands can be grouped into two different types according to their bulk major elements chemistry: (i) OI1-u depleted composition (xenoliths from La Palma, Lanzarote and Tenerife HEXO, Table 1); (ii) OI2-u depleted & metasomatized composition (Tenerife HLCO, Table 1). Geophysical parameters are affected to a different extent by changes in the bulk mantle composition for given temperature pressure conditions. Fig. 3 illustrates the impact on the calculated mantle densities and elastic seismic velocities of varying the bulk composition for a reference 100-km-thick lithospheric column. Depleted compositions yield similar mantle properties, showing an increasing (decreasing) trend for the density (seismic velocities) going from East – Lanzarote – to West – La Palma. This trend mimics the degree of compositional depletion (i.e., increasing Mg#) from La Palma to Lanzarote (Table 1). Mantle composition from Tenerife HEXO xenoliths lies approximately in between La Palma and Lanzarote xenolith compositions, with associated densities and seismic velocities in the middle of the range for OI1-u compositions. OI2-u composition (Tenerife HLCO) is associated with comparatively high densities similar to those of an

average ocean peridotite and Azrou Middle Atlas xenoliths compositions, although the spinel-garnet phase transition (density/velocity discontinuity) at around 40 km depth is missing for Tenerife HLCO composition (Fig. 3). 6.2. Purely lithospheric models Two end-member models can be defined in terms of the mantle composition assumed in domain 2 (Fig. 1). Models C1 and C2 are defined by OI1-u (Lanzarote) and OI2-u (Tenerife HLCO) mantle compositions in the Canary Islands respectively (Table 1). The assumed OI1-u composition is specifically based on xenoliths from Lanzarote in C1. The rationale is that among the depleted, preocean-island compositions listed in Table 1, the density associated with the xenoliths from Lanzarote is the lowest (Fig. 4). Therefore, C1 and C2 represent end-member lithospheric density distributions beneath the Canary Islands. The sub-lithosphere in C1 and C2 is compositionally homogeneous (PUM, Table 1) and without major deep thermal anomalies (i.e., only small temperature variations given by the lateral variations in the lithospheric thickness). The lithospheric thickness in model C1 and C2 required to fit the measured geophysical/geodatic observables (Fig. 2), given the assumed crustal structure (see Appendix A), differs in 20–40 km, with the maximum discrepancies being located beneath the central and eastern islands (Fig. 4). Model C2 shows a moderate lithospheric thinning beneath the Canary archipelago (LAB depth of 90–100 km) compared to the 110–130-km-thick lithosphere in the Jurassic oceanic and continental Moroccan domains (Fig. 4). In contrast, the lithosphere in the archipelago in model C1 is as thick as in the continental area and slightly thicker than in the oceanic domain. Both C1 and C2 models allow fitting reasonably well the long-wavelength part of most of the constraining observables, excepting elevation calculated according to local isostasy (Fig. 5, Table 2). The elevation misfits are discussed in detail in Section 6.4 in the light of a flexural isostatic analysis. Further insight can be gained by comparing models C1 and C2 against seismic tomography. The multimode inversion of surfaceand S-wave forms in Legendre et al. (2012) fits one of the objectives of our study: compare the amplitudes of synthetic and

J. Fullea et al. / Earth and Planetary Science Letters 409 (2015) 71–88

77

Fig. 4. (A) to (D) depth of the thermal lithosphere–asthenosphere boundary in models C1, C2, C3 and C4 respectively (see the text for further details).

tomography seismic velocity anomalies to infer the temperature distribution in the sub-lithospheric mantle. In contrast to bodywave tomography models, for which the amplitude of velocity anomalies is generally not well resolved, surface-wave inversions provide robust estimates of the absolute shear-wave values which can be thoroughly interpreted in terms of the mantle’s thermochemical structure (e.g., Connolly and Kerrick, 2002; Afonso et al., 2008). To make an appropriate comparison, elastic mantle velocities computed from our thermochemical models have to be corrected for attenuation effects which are intrinsically included in seismic data. Here, anelasticity is computed using a pressure- and temperature-dependent formula (e.g., Minster and Anderson, 1981; Karato, 1993; Afonso et al., 2005; Fullea et al., 2012) for the relevant reference period (50 s for the surface-wave tomography study of Legendre et al., 2012) and a reference grain size of 10 mm. To include the intrinsic vertical smearing in surfacewave tomography models, we consider vertical averages of the seismic models over broad depth ranges roughly corresponding to the lithosphere (45–130 km) and sub-lithosphere (130–295 km). Fig. 6 shows a vertical average of the synthetic seismic velocity anomalies predicted by models C1 and C2 and the model of Legendre et al. (2012) at lithospheric and sub-lithospheric depth ranges for a comparison. At lithospheric depths model C1 shows a positive velocity anomaly, southwestwards of the Canary archipelago, which is partially coincident with a negative anomaly centered (MINECO) according to Legendre et al. (2012). Synthetic velocity anomalies in model C2 are more similar to the seismic tomography model of Legendre et al. (2012) in that the archipelago is underlain by a negative anomaly, although this anomaly has a clear NE–SW trend along the North Atlantic African

margin that is missing in Legendre et al. (2012). Both model C1 and C2 fail to reproduce the observed sub-lithospheric seismic anomaly over the Canary Islands. The reason for it is that models C1 and C2 are nearly homogeneous in the sub-lithosphere, i.e. without major deep thermal and/or compositional anomalies. 6.3. Sub-lithospheric anomaly The existence of a low velocity zone in sub-lithospheric mantle below the Canary archipelago, as imaged by seismic tomography (e.g., Hoernle et al., 1995; Legendre et al., 2012; Bonnin et al., 2014), can be simulated by introducing an ad hoc thermal and/or compositional anomaly under the lithosphere in the thermochemical models (C1 or C2). In the context of our modelling, synthetic mantle seismic velocities in the models can be reduced by either mantle heating or through compositional changes (e.g., diminishing Mg#, see Fig. 3). Considering that the sub-lithospheric composition assumed in our models (PUM, Table 1) is already pristine and, in this sense, very fertile (i.e., low velocities), any compositional variation would go into the direction of getting a more refractory mantle which would be potentially faster (and less dense) than the ambient PUM (cf. Fig. 3). Hence, compositional changes are unlikely the cause of the negative/low sub-lithospheric velocity anomaly imaged beneath the Canary Islands (Fig. 6). On a broader context, small amounts of partial melt could be, along with a high temperature anomaly, the cause of a low velocity anomaly in the mantle, although analyzing melt effects is out of the scope of this paper. The mantle potential temperature, T pot , in our models corresponds to around 1335 ◦ C in the sub-lithosphere. This is well

78

J. Fullea et al. / Earth and Planetary Science Letters 409 (2015) 71–88

Fig. 5. Residuals (synthetic − measured) for the lithospheric model C1 (see Table 2). (A) Geoid anomaly. (B) Free air anomaly. (C) Bouguer anomaly. (D) Elevation (local isostasy). (E1–E6) Gravity gradients computed at 255 km above the sea level (GOCE mission perigee height). The sub-indexes of U denote second derivatives of the Earth’s potential (i.e., gravity gradient) along the west–east, south–north and vertical directions (x, y, z respectively).

J. Fullea et al. / Earth and Planetary Science Letters 409 (2015) 71–88

Table 2 Statistics of the residuals (synthetic-measured) for the constraining geophysical and geodetic observables used in this study (Fig. 2). Std. dev. residuals

Model C1a

Model C2a

Model C3a

Model C4a

Topo (m) FA (mGal) Bouguer (mGal) Geoid (m) Uzz (mE) Uxx (mE) Uyy (mE) Uzx (mE) Uzy (mE) Uxy (mE)

382 14.01 14 1.01 74.6 40.2 68.5 74.8 96.6 71.3

367.3 12.04 11.77 0.77 61.4 80 69.9 101.04 80.1 87.3

374.5 12.2 11.8 1 .9 70.7 81.2 69.2 138.2 105.3 67.2

444.1 13.03 13.2 1.43 78.5 79.6 90.1 182.1 100.7 61.2

a For details on the mantle compositions in models C1–C4 see Table 1 and the text. Details on the source of measured data are given in Section 4 of the text. In all models, the gravity gradients are computed at the satellite altitude (255 km). The thermal characteristics of models C1–C4 are summarized in Table 3.

within the expected range of T pot = (1280–1400) ◦ C in ambient mantle around the world excluding mantle plumes (e.g., Herzberg et al., 2007; Herzberg and Asimow, 2008). Other well studied hot spots like Hawaii or Iceland are characterized by T pot of 1615–1790 ◦ C based on olivine-liquid equilibria, a relatively high value compared to standard mid-ocean ridges (1450 ◦ C–1475 ◦ C) (Putirka, 2005; Putirka et al., 2007). In the modelling approach followed in this study (LitMod-based) isostatic topography (either local or flexural) is assumed to be static, i.e. no dynamic contributions from the convecting mantle are implicitly included (e.g., Afonso et al., 2008; Fullea et al., 2009). Therefore, any thermal anomaly in the sub-lithospheric mantle in our models will not have, through the corresponding density anomaly, a direct impact on the calculated static elevation. Yet sub-lithospheric density anomalies associated with thermal anomalies introduced ad hoc in the models do have an impact in other density-dependent observables (gravity and geoid anomalies, gravity gradients). That leads to an indirect effect of sub-lithospheric thermal anomalies on the predicted static elevation: if a purely lithospheric model (e.g., C1 or C2) able to explain the measured data (excepting elevation, see Fig. 5 data and Table 2) is modified to include a temperature anomaly in the sub-lithosphere, then the predicted data will no longer match the observed values and further modifications to the lithospheric structure will be required to counterbalance the new sub-lithospheric density distribution. Those changes required in the lithospheric density distribution (e.g., thickness of the lithosphere) will, in turn, affect the calculated static elevation (i.e. buoyancy): hence the indirect effect over elevation of sub-lithospheric thermal/density anomalies. This indirect effect explains why synthetic (static) topography differs among models with the same crustal structure and mantle composition but different sub-lithospheric temperature/density distributions. Model C3 represents a modification of model C2 (i.e. OI2-u composition in the Canarian domain, Fig. 1) where a sublithospheric thermal anomaly of +100 ◦ C with respect to the ambient mantle has been introduced ad hoc from the base of the lithosphere down to the base of the model at 400 km depth in parts of lithospheric domain 2 (Fig. 7). Variations in the lithospheric thickness to compensate for the deep temperature/density anomaly are allowed so as the measured observables are matched (similar residuals for models C1, C2 and C3, Table 2). T pot in model C3 is 1435 ◦ C compared to T pot = 1335 ◦ C for purely lithospheric models C1 and C2 (i.e. no sub-lithospheric anomalies included). The lithospheric thickness in model C3 is moderately different from those of models C1 and C2 (differences <30 km) (Fig. 4). In contrast, synthetic sub-lithospheric velocity anomalies in C3 show a negative anomaly similar in amplitude (−(2–2.5)%) to that in the model of Legendre et al. (2012) (Fig. 6).

79

A useful exercise is to re-compute model C3 with a thermal sub-lithospheric anomaly of +300 ◦ C, instead of +100 ◦ C, with respect to the ambient mantle. In model C4 (T pot = 1635 ◦ C) the corresponding sub-lithospheric low density anomaly (≈−30 kg/m3 ) is larger than in the case of model C3 (Fig. 7), and, hence, larger modifications of the lithospheric structure are required to fit most of the measured data (Table 2). Assuming OI1-u or OI2-u compositions as the average bulk composition in the Canary Islands domain (Fig. 1, Table 1), the lithosphere under the archipelago needs to be modelled thicker than in the neighbor oceanic and Moroccan continental domains (>130 km depth; Fig. 4). Furthermore, the predicted synthetic negative velocity anomaly beneath the Canarian domain in the sub-lithosphere for model C4 is significantly lower (≈−8%) than those of model C3 and estimated by Legendre et al. (2012) based on surface-wave data (≈−2.5%). The residuals (calculated − measured observables) for model C4 are similar to those of models C1, C2 and C3 for gravity-fieldderived geophysical and geodetic observables (geoid and gravity anomalies, and gravity gradient) (Table 2). The variations in the residuals for elevation differ considerably among the thermochemical models presented here and are discussed in the following section. 6.4. Flexural isostasy All the thermochemical models presented in this work (C1, C2, C3 and C4) show a static residual elevation (synthetic-observed), computed assuming local isostasy (see Section 5.4), that is clearly negative beneath the Canaries, i.e. the synthetic elevation is lower than the actual topography of the islands (Fig. 5D). This misfit is to be expected given the relative short wavelength and high topography of the islands. Hence, a further flexural isostatic calculation, where part of the topographic load is supported by the strength of the lithosphere, is required for an appropriate analysis of the isostatic balance. In our approach, the flexural load is given by the departure of the synthetic, locally isostatic elevation from the measured topography/bathymetry (see Section 5.4). Fig. 8 shows the results of the isostatic flexural balance for the thermochemical models presented here. The case where T e = 0 km (Fig. 8, A1–A3) corresponds to a local isostatic balance under the assumption that shear stresses across the vertical boundaries of every (ideal) lithospheric column are negligible (i.e., free vertical movement). In models C1 and C2 the residual locally isostatic elevation (synthetic-observed) in the Canarian lithospheric compositional domain is positive (except over the islands) and around +500 m on average. The residual elevation for T e = 0 km is considerably lower in most parts of the domain for model C3, and is negative (≈−500 m on average) in model C4. For a mildly weak lithosphere with T e = 20 km (e.g., thermally rejuvenated old oceanic lithosphere) the residual flexural isostatic elevation is still negative for the islands in all models, whereas the surrounding areas (swell and margin) exhibit negative (≈−400 in C4, north of the archipelago) to positive (100–600 m in C1, and C3) values (Fig. 8, B1–B3). For an intermediate lithospheric rigidity (T e = 35 km) the situation is similar to the case in which T e = 20 km, with somewhat smaller amplitudes in the positive and negative residual flexural elevations (Fig. 8, C1–C3). If high T e values are considered (i.e., T e = 85 km, corresponding to mechanically strong lithosphere), the residual topography is almost zero over the archipelago for models C1 and C2, with a smooth positive residual south of the islands (≈150 m). In contrast, there is a long wavelength negative residual flexural elevation along the islands, generally more pronounced over the western than over the eastern islands, in models C3 (−(150–400) m) and C4 (−(600–800) m) (Fig. 8, D1–D3).

80

J. Fullea et al. / Earth and Planetary Science Letters 409 (2015) 71–88

Fig. 6. Seismic tomography models. (A) and (B) Multimode inversion of surface- and S-wave forms in Europe from Legendre et al. (2012). (A) V s anomalies averaged between depths of 45–130 km. (B) V s anomalies averaged between depths of 130–295 km. (A1–A4) synthetic S-wave velocity anomaly averaged between depths of 45–130 km for models C1, C2, C3 and C4 respectively. (B1–B4) synthetic S-wave velocity anomaly averaged between depths of 130–295 km for models C1, C2, C3 and C4 respectively. The corresponding average velocity at each depth has been subtracted to compute the synthetic anomalies in A1–A4 and B1–B4. Anelasticity effects are included in the synthetic seismic velocity anomalies considering the dominant frequencies of the seismic tomography model of Legendre et al. (2012) (∼50 s) and an average grain size of 1 cm (see text for further details).

J. Fullea et al. / Earth and Planetary Science Letters 409 (2015) 71–88

81

Fig. 7. Upper mantle S-wave velocity and density distribution in models C2, C3 and C4. (A) Elevation map with lithospheric mantle compositional domains as in Fig. 1. (B), (C) and (D) synthetic S-wave velocity anomaly averaged between depths of 130–295 km for models C2, C3 and C4 respectively. The solid black line in (A–D) shows the location of the V s and density profiles displayed below. (E1), (F1), (G1) Synthetic S-wave profiles for models C2, C3 and C4 respectively. (E2), (F2), (G2) Density profiles for models C2, C3 and C4 respectively. The crustal layers are shown on top of the lithosphere (see Appendix A for a detailed description of the crustal model). The black solid lines are the mantle compositional boundaries (see A in this figure for location) and the thermal LAB (see Fig. 4). In E1 and E2 the LAB of model C1 is shown for a comparison (black dashed line).

7. Discussion and conclusions In this work, the thermal and compositional structure of the lithosphere and uppermost mantle beneath the Canary archipelago has been assessed by means of an integrated geophysical– petrological approach able to self-consistently combine various data sets (Afonso et al., 2008; Fullea et al., 2009). A number of alternative thermochemical models (C1–C4) have been presented here based on geophysical/geodetic (i.e., elevation, gravity and geoid anomalies, gravity gradients, seismic) and mantle xenolith constraining data (Table 1). These models can be grouped into two broad categories: C1 and C2 are purely lithospheric models whereas C2 and C3 include sub-lithospheric thermal anomalies

of variable magnitude reflecting the commonly accepted range of hot spot high anomalous temperatures with respect to ambient mantle, +(100–300) ◦ C (e.g., Herzberg and Asimow, 2008; White, 2010) (Table 3). 7.1. Mantle composition beneath the Canary archipelago The major element chemistry of the mantle xenoliths erupted in the Canary Islands (see Section 3.2) reflects a remarkably refractory mantle as a result of intense melt extraction (up to 25–30%; Neumann et al., 2004). The ultra-depleted nature of the mantle xenoliths in the Archipelago seems to be a common feature in most ocean islands around the world (Simon et al., 2008). Models

82

J. Fullea et al. / Earth and Planetary Science Letters 409 (2015) 71–88

Fig. 8. Flexural isostatic residual elevation (synthetic-measured) for different elastic thicknesses, T e , in models C1, C3 and C4 (see the text for further details). (A1–A3) Residuals assuming T e = 0 km for models C1, C3 and C4 respectively. (B1–B3) Residuals assuming T e = 20 km for models C1, C3 and C4 respectively. (C1–C3) Residuals assuming T e = 35 km for models C1, C3 and C4 respectively. (D1–D3) Residuals assuming T e = 85 km for models C1, C3 and C4 respectively. We note that elevation residuals for purely lithospheric models C1 and C2 are very similar and therefore only results for model C1 are shown. Table 3 Summary of the thermochemical models presented in this study.

Mantle compositiona T pot (◦ C)b Lithospheric thickness (km)c a b c

Model C1

Model C2

Model C3

Model C4

OI1-u (Lanzarote) 1335 110–130

OI2-u (HLCO Tenerife) 1335 90–100

OI2-u (HLCO Tenerife) 1435 100–110

OI2-u (HLCO Tenerife) 1635 130–140

Mantle composition in the Canarian lithospheric domain (see Fig. 1); for details on the bulk amount of major oxides see Table 1 and the text. Mantle potential temperature. Average LAB depth in the Canary Islands, see Fig. 4 and Fig. 7.

C1 and C2 are compositional and density end-member lithospheric models: C1 is defined by OI1-u (Lanzarote, Table 1) depleted composition, likely representative of the upper mantle composition prior to the onset of Canary Islands magmatism, whereas man-

tle composition in C2 is OI2-u (Tenerife HLCO, Table 1) hence reflecting subsequent metasomatism (silicic carbonatite melt). As pointed out by Simon et al. (2008), OI2-u composition similarities with other relatively fertile compositions (e.g., average ocean

J. Fullea et al. / Earth and Planetary Science Letters 409 (2015) 71–88

peridotite and Middle Atlas xenoliths, Table 1) does not imply the former has to be associated with pristine mantle or mantle that has undergone lower degrees of partial melting than OI1-u harzburgites: OI2-u is to be interpreted as the result of the interaction between OI1-u protoliths and melts. Yet from a geophysical perspective C2 can be regarded as a compositionally homogeneous model, i.e. a unique mantle domain with a relatively fertile or depleted and metasomatized composition, due to the similar density associated with compositions OI2-u, average ocean peridotite and Middle Atlas (Fig. 3). Differences in the thermal thickness of the lithosphere between models C1 and C2 are moderated (<40 km, Fig. 4). Interestingly, in neither model C1 nor C2 the lithosphere beneath the Canaries is particularly thinned (90–130 km, Fig. 4) with respect to the surrounding continental and oceanic domains (Table 3). Similar values are obtained by Ranero et al. (1995) in contrast to the weak and thermally rejuvenated lithosphere inferred by Watts (1994) and Watts et al. (1997) based on gravity and seismic T e estimates. 7.2. Mantle potential temperature Absolute values for potential temperatures in the mantle beneath hot spots can differ by >100 ◦ C according to different authors (cf. Herzberg, 2011). Putirka et al. (2007) have estimated potential temperatures of 1616 ◦ C, 1722 ◦ C and 1454 ± 81 ◦ C for Iceland, Hawaii and average mid-oceanic ridge respectively, based on olivine-liquid thermometers. In contrast, Herzberg and Asimow (2008) reduce these estimates to 1450 ± 10 ◦ C, 1560 ± 40 ◦ C and 1350 ◦ C for Iceland, Hawaii and an average mid-oceanic ridge respectively, on the basis on melting modelling of observed lava compositions. Yet there is a consensus in that ocean island hot spots cannot be explained solely by volatile enrichment or source fertility: there must be a thermal anomaly with respect to ambient mantle, and this anomaly is likely in the range of +(100–300) ◦ C (e.g., Herzberg and Asimow, 2008; White, 2010). According to Putirka (2008), the potential temperatures for ambient mantle and in the Canary Islands are 1396 ◦ C and 1561 ◦ C respectively. Models C3 and C4 include sub-lithospheric thermal anomalies of +100 ◦ C and +300 ◦ C with respect to ambient mantle (here assumed to be 1335 ◦ C) and illustrate therefore the end member assumed mantle temperature anomalies in hot spots (Table 3). Lithospheric thickness in model C3 lies somewhere in between the values for models C1 and C2, showing differences of <30 km (Fig. 4). Model C4 (T pot = 1635 ◦ C) can be ruled out as a candidate to represent the present day thermal structure beneath the Canary archipelago: (i) the amplitude of sub-lithospheric predicted V s anomalies in the depth range (130–295) km is too high (≈−8%) with respect to surface-wave tomography values (≈−2.5%); and (ii) the amplitude of the residual flexural isostatic elevation in the Canaries is too large (i.e., negative buoyancy) regardless the assumed effective elastic thickness value (Fig. 8). Models C1 and C2 (T pot = 1335 ◦ C) do not show any significant sub-lithospheric seismic anomaly and are at odds with tomographic images (Fig. 6). Model C3 is able to explain both observed seismic anomalies (Fig. 6) and measured geophysical/geodetic data (gravity and geoid anomalies, gravity gradients, Table 2). The potential temperature in C3 (1435 ◦ C) coincides with lava-composition-derived estimates of 1420–1480 ◦ C, within the range of expected values for ocean islands in the Atlantic (Herzberg and Asimow, 2008), but is >100 ◦ C colder than T pot = 1561 ◦ C derived by Putirka (2008) for the Canary archipelago. According to Putirka et al. (2007) a potential temperature of 1435 ◦ C in the Canary archipelago, as in model C3, would fall within the typical range for ambient mid-ocean ridges, instead of ocean islands. Our results favor an ambient mantle temperature of 1335 ◦ C as in Herzberg et al. (2007) and Herzberg and Asimow (2008) rather than 1400–1475 ◦ C (Putirka et al., 2007;

83

Putirka, 2005, 2008). In this work the low velocity anomaly in the mantle beneath the Canary Islands is fully explained by a high temperature anomaly as first order approximation. However, small packets of melt distributed in pore spaces cannot be ruled out as an additional cause of low shear-wave values in the mantle. For the sake of simplicity, here the Canary thermal anomaly is modelled as a sharp temperature contrast in the mantle (Fig. 7); in the real Earth, though, the core of the thermal anomaly would be hotter than its periphery, where the presence of volatiles (e.g., CO2 and H2 O) would likely control mantle melting processes. Therefore, the periphery of the low velocity zone modelled here could be caused, to some extent, by low degree CO2 -rich melts rather than only high temperatures with respect to ambient mantle. 7.3. Strength of the lithosphere A number of studies have discussed the need to estimate the flexural response of the Canary Islands topographic load; in addition, the amplitude of the Canary swell is arguably a key part of the discussion on the origin of the islands (see references in Section 3.3). The range of effective elastic thickness (T e ) values estimated for the archipelago varies from 20 km (Watts, 1994; Watts et al., 1997) to 50 km (Filmer and McNutt, 1989). Here we consider, for the whole study area, a range of constant and isotropic T e values spanning from 0 km (i.e., local isostasy) to 85 km (i.e., mechanically strong lithosphere). It is interesting to note that for a weak lithosphere (T e ≤ 20 km) the residual isostatic flexural elevation (calculated − observed) is still negative (−(400–1400) m) over the islands in all models (C1–C4) (Fig. 8, B1–B3). The density distributions in models C1–C4 are able to reproduce the long wavelength component of the other observables (Fig. 5, Table 2) but, at the same time, predict a negative buoyancy beneath the islands and a positively buoyant area in the surroundings (swell and margin) if low T e values (<20 km) are assumed. The former implies that the hypothesis of a thermally rejuvenated (weakened) old oceanic lithosphere is unlikely based on the results presented here. In contrast, our results suggest that the lithosphere should be mechanically strong, in line with the range of thermal lithospheric thicknesses obtained in this study (Fig. 4). The residual flexural elevation for models including a sublithospheric thermal anomaly (C3 and C4) is always negative (Fig. 8). A negative value in the residual elevation (under either local or flexural isostasy) means that the static buoyancy of the lithosphere in the model is lower than expected. For a strong lithosphere (T e of 85 km), consistent with the modelled lithospheric thicknesses (Fig. 4), this residual attains unacceptably high values of −(600–800) m in model C4. In contrast, the residuals for model C3 are relatively moderate (Fig. 8). We deem this value of 150–400 m a required dynamic contribution to the static predicted elevation in C3 if the present-day observed topography in the Canary Islands and associated swell is to be reproduced. Such dynamic contribution would be linked to the sub-lithospheric thermal anomaly of +100 ◦ C (and ≈−10 kg/m3 ) with respect to the ambient mantle (Fig. 7). The mechanically strong lithosphere (i.e., high viscosity) suggested by our models could explain the lack of a broad topographic swell as observed in other hot spots (e.g., Hawaii), in spite of the dynamic elevation related to the sub-lithospheric thermal anomaly (i.e. mantle upwelling). Models C1 and C2 do not include any significant sub-lithospheric thermal/density anomaly implying that no dynamic topography is expected. The relatively uniform lithospheric thickness below the Canary Islands and adjacent continental domain questions the existence of small-scale edge driven convection interacting with a mantle plume, as proposed by some authors to explain the irregular space-temporal volcanic pattern in the islands (e.g., Geldmacher et al., 2005).

84

J. Fullea et al. / Earth and Planetary Science Letters 409 (2015) 71–88

7.4. Conclusions – Compositional differences in the mantle based on ultrarefractory xenoliths from the Canary Islands translate into moderate thermal LAB variations beneath the archipelago (<40 km). The thermal lithosphere beneath the islands is not thin compared to the surrounding oceanic and continental domains (110 ± 20 km). – A model (C3) including a thermal sub-lithospheric anomaly of +100 ◦ C (T pot = 1435 ◦ C) beneath the Canaries is able to explain both observed seismic velocity anomalies (in the sublithosphere) and measured geophysical and geodetic data (Table 3). Mantle potential temperatures as high as 1635 ◦ C would predict: (i) V s anomalies in the sub-lithosphere too negative with respect to surface-wave tomography, and (ii) a too negatively buoyant lithosphere regardless the assumed T e value. – Models without a positive temperature anomaly in the sublithosphere (C1 and C2, characterized by mantle T pot = 1335 ◦ C, see Table 3) fail to reproduce the observed sub-lithospheric seismic anomaly over the Canary Islands. – Based on a flexural analysis of the static elevation residuals the hypothesis of a thermally rejuvenated (weakened) old oceanic lithosphere (T e ≤ 20 km) is unlikely. Our models suggest that the lithosphere is mechanically strong, in line with the range of thermal lithospheric thicknesses obtained in this study. This suggests that significant thermal erosion at the base of the lithosphere did not take place coeval to the most recent magmatic activity. – A sub-lithospheric anomaly of +100 ◦ C (model C3) requires a contribution of 150–400 m of dynamic topography to the static predicted elevation to match the present-day observed topography in the Canary Islands and associated swell. Acknowledgements Javier F was supported by the JAE-DOC programme (CSICSpain) cofounded by ESF, and by Spanish Ministry of Economy and Competitiveness in Gran Canaria grants CGL2009-13103 and CGL2012-37222. José F and AGC have been supported by MINECO research project AYA2010-17448. This is a contribution for the Moncloa Campus of International Excellence. We want to thank the editor (J.P. Brodholt), C. Herzberg and one anonymous reviewer for their constructive comments and suggestions which helped improving the manuscript. Appendix A. Crustal model An initial crustal model is derived from 1D inversion of geoid anomaly and elevation data (Fullea et al., 2006, 2007). The method considers a simple two-layer lithospheric model (i.e. crust and lithospheric mantle) in which crustal density varies linearly with depth whereas lithospheric mantle density is temperaturedependent only (i.e., neither compositional nor pressure effects are included). In spite of its limitations (i.e., 1D approach, no compositional effects), the inversion of geoid anomaly and elevation data is a useful tool to define a starting, reference crustal model. This initial model has been subsequently modified to build a more comprehensive model integrating available seismic constraints in the Canary archipelago and neighboring Atlantic African passive margin as summarized below. Available seismic information covers the central and eastern islands and the Atlantic Moroccan passive margin (see Fig. 1). The oceanic crust around the archipelago exhibits a normal velocity structure with a thickness of 5–7 km (e.g., Banda et al., 1981). Beneath the central islands (Gran Canaria and Tenerife) the

crust comprises, according to wide-angle seismic refraction studies, a 10–12 km-thick upper volcanic edifice (V p = 5.5–6 km/s) and a lower crust (6.6–7.3 km/s), which is 2–3 km thick offshore and thickens beneath the islands to around 6 km (Ye et al., 1999; Dañobeitia and Canales, 2000; Watts et al., 1997). In Gran Canaria the lower volcanic edifice consists of a 4–6 km thick low velocity zone (V p = 5.5 km/s), interpreted as an old Miocene syenitic feldspar-rich core (Ye et al., 1999). In the eastern islands (Lanzarote and Fuerteventura) the uppermost layer (4–5 km) shows an average V p of 4 km/s, whereas upper crust varies in thickness (4–6 km in Lanzarote and 10–12 km in Fuerteventura) and P-wave velocity (6 km/s in Lanzarote and 6.5–6.7 km/s in Fuerteventura) (Dañobeitia and Canales, 2000 and references therein). Analysis of P- and S-wave receiver functions (RF) generated on seismic interfaces in the archipelago reveals crustal discontinuities at depths of 18–22 km in Lanzarote, 23.5 km in Gran Canaria, 14–18 km in Tenerife and 11–14 km beneath La Palma and el Hierro (Lodge et al., 2012; Martinez-Arevalo et al., 2013). These RF’s estimates are compatible with the Moho (V p = 7.4 km/s) derived from seismic experiments in Tenerife (Dañobeitia and Canales, 2000) and coincide with the base of a magmatic underplating level (V p = 7.4–8 km/s, see below) in Gran Canaria (Ye et al., 1999; Krastel and Schmincke, 2002) and Lanzarote (e.g., Dañobeitia and Canales, 2000). The Atlantic Moroccan passive margin north of the Canary archipelago (latitude 33◦ N) is characterized by a progressive transition from old oceanic crust to a 34–35 km-thick continental crust consisting on three layers: upper (V p = 5.2–6.4 km/s), middle (V p = 6.6–6.8 km/s) and lower crust (V p = 6.8–7 km/s) (SISMAR profile, Contrucci et al., 2004). South of the Canary Islands (latitude 22◦ N) the Atlantic Moroccan continental margin lacks the lower crustal layer present in its northern part (Moho depth of about 27 km) and its seismically described by 27-km-thick crust with an upper (V p = 6 − 6.2 km/s) and a middle (V p = 6.5–6.8 km/s) layer (Dakhla profile, Klingelhoefer et al., 2009). A recent teleseismic P- and S-wave RF’s study in the western Moroccan Atlas shows two stations, MM13 and MM14, close to the Atlantic margin (north of our study region) with Moho depth estimates of 24.65 ± 6.7 km and 28.05 ± 6.7 km respectively (Spieker et al., 2014). Underlying the lower crust Ye et al. (1999) interpret an 8–10 km-thick layer of magmatic underplating in Gran Canaria, defined by a range of P-wave velocities (7.4–8 km/s) higher than those of typical lower crustal material, but lower that the average Jurassic oceanic uppermost mantle velocities in the neighborhood of the islands (8 km/s). A reinterpretation of the seismic sections in Banda et al. (1981) and profile CD82-P12 (Tenerife) led Dañobeitia and Canales (2000) to infer the existence of a layer with similar anomalous velocity (7.4–8 km/s) beneath Tenerife (11 km-thick) and the eastern islands (7–12 km-thick). Lodge et al. (2012) also found evidence for magmatic underplating beneath Lanzarote (8 km), La Palma (4 km) and Tenerife (2–5 km) based on RF’s analysis. Further petrological support for the magmatic underplating as a cumulate of basic magma impinging at the base of the crust in Gran Canaria is provided by the evidence of a reservoir at around 14 km depth, where fractionation of felsic magma took place (Freundt and Schmincke, 1995). The radial extension of the magmatic material added at the base of the crust is rather limited, at least in central islands, where the underplating ends abruptly around 20 km away from the shoreline (Ye et al., 1999; Dañobeitia and Canales, 2000). In view of the structural complexity of the area we define a crustal model able to represent three broad tectonic crustal domains based on the available seismic information: oceanic, volcanic island and continental crust. The oceanic crust includes sediments and a thin igneous layer with V p values of 5.2–7.3 km/s with an average density of 2940 kg/m3 . The Canary Islands are defined by a volcanic edifice/upper crust (5.2–6.4/5.5–6 km/s), a lower crust

J. Fullea et al. / Earth and Planetary Science Letters 409 (2015) 71–88

85

Fig. 9. Crustal model. (A1) Elevation with the vertical axis positive downwards. (A2–A6) Depth to the base of the different crustal layers. (A2) layer 1: sediments, (A3) layer 2: upper crust/volcanic edifice, (A4) layer 3: middle/oceanic crust, (A5) layer 4: lower crust and (A6) layer 5: magmatic underplating (Table 4). Each crustal layer is bounded by its base and the base of the layer above. Table 4 Geophysical properties of the different crustal bodies used in the 3D models (C1, C2, C3 and C4). See Figs. 9 and 10 for further details on the geometry of the crustal layers. Layer

Density (kg/m3 )

Vp (km/s)

Heat production (μW/m3 )

Thermal conductivity (W/m K)

1) 2) 3) 4) 5)

2300 2670 2940 3020 3200

<4.5 5.2–6.4/5.5–6 6.6–6.8/5.2–7.3 6.6–7.3 7.4–8

1.5 1 0.5 0.2 0.1

2.5 2.5 2.5 2.1 2.1

Sediments Upper crust/volcanic edifice Middle/oceanic crust Lower crust Magmatic underplating

(6.6–7.3 km/s) and a layer of magmatic underplating (7.4–8 km/s). Continental African crust comprises sediments, upper, middle and lower crust. Therefore, up to five layers are defined in our crustal model: (1) sediments, (2) upper crust/volcanic edifice, (3) mid-

dle/oceanic crust, (4) lower crust and (5) magmatic underplating (Fig. 9, Fig. 10, Table 4). This number of layers represents a compromise between the complexity of the study area and the need to reduce the number of layers to build a manageable model.

86

J. Fullea et al. / Earth and Planetary Science Letters 409 (2015) 71–88

Fig. 10. Crustal model. (A) Moho depth map. The black dotted lines are the locations of the crustal profiles P1–P4 displayed below. (B) Thickness map of the magmatic underplating. P1–P4 are representative cross-sections showing the crustal structure.

Two crustal layers are grouped according to their similar average V p values: layer 2 plays de role of either an upper crust or a volcanic edifice; layer 3 describes the middle continental and oceanic crusts. Densities of the crustal layers are computed based on V p -density relationships (Christensen and Mooney, 1995) (Table 4). An appropriate layering of the available crustal units allows us to define the density stratification in different crustal tectonic domains: oceanic domain → sediments + oceanic crust, continental crust domain → sediments + upper + middle + lower crust, volcanic islands domain → sediments + volcanic edifice + lower crust + magmatic underplating (Fig. 10). Although seismic constraints suggesting the presence of magmatic underplating are only available for the central and eastern islands and La Palma (e.g., Dañobeitia and Canales, 2000), we have extended it to the other islands (La Gomera and El Hierro) and north of Lanzarote based on our modelling results (Figs. 9 and 10). References Ablay, G.J., Kearey, P., 2000. Gravity constraints on the structure and volcanic evolution of Tenerife, Canary Islands. J. Geophys. Res. 105, 5783–5796. Afonso, J.C., Zlotnik, S., 2011. The subductability of the continental lithosphere: the before and after story. In: Brown, D., Ryan, P.D. (Eds.), Arc-Continent Collision. In: Frontiers in Earth Sciences. Springer, pp. 53–86. Afonso, J.C., Fernàndez, M., Ranalli, G., Griffin, W.L., Connolly, J.A.D., 2008. Integrated geophysical–petrological modeling of the lithosphere and sublithospheric upper mantle: methodology and applications. Geochem. Geophys. Geosyst. 9, Q05008. http://dx.doi.org/10.1029/2007GC001834. Afonso, J.C., Ranalli, G., Fernàndez, M., 2005. Thermal expansivity and elastic properties of the lithospheric mantle: results from mineral physics of composites. Phys. Earth Planet. Inter. 149, 279–306.

Ancochea, E., Barrera, J., Bellido, F., Benito, R., Brändle, J., Cebriá, J., Coello, J., Cubas, C., De La Nuez, J., Doblas, M., 2004. Canarias y el vulcanismo neógeno peninsular. In: Vera, J.A. (Ed.), Geología de España. SGE-IGME, Madrid, pp. 635–682. Anguita, F., Garcia Cacho, L., Colombo, F., González Camacho, A., Vieira, R., 1991. Roque Nublo caldera: a new stratocone caldera in Gran Canaria, Canary Islands. J. Volcanol. Geotherm. Res. 47, 45–63. Anguita, F., Hernán, F., 2000. The Canary Islands origin: a unifying model. J. Volcanol. Geotherm. Res. 103, 1–26. Araña, V., Camacho, A., Garcia, A., Montesinos, F., Blanco, I., Vieira, R., Felpeto, A., 2000. Internal structure of Tenerife (Canary Islands) based on gravity, aeromagnetic and volcanological data. J. Volcanol. Geotherm. Res. 103, 43–64. Banda, E., Dañobeitia, J.J., Suriñach, E., Ansorge, J., 1981. Features of crustal structure under the Canary Islands. Earth Planet. Sci. Lett. 55, 11–24. Bonnin, M., Nolet, G., Villaseñor, A., Gallart, J., Thomas, C., 2014. Multiple-frequency tomography of the upper mantle beneath the African/Iberian collision zone. Geophys. J. Int. 198 (3), 1458–1473. Bowin, C., 2000. Mass anomalies and the structure of the Earth. Phys. Chem. Earth 25, 343–353. Camacho, A., Fernández, J., González, P., Rundle, J., Prieto, J.F., Arjona, A., 2009. Structural results for La Palma island using 3-D gravity inversion. J. Geophys. Res. 114, B05411. http://dx.doi.org/10.1029/2008JB005628. Camacho, A., Fernández, J., Gottsmann, J., 2011. A new gravity inversion method for multiple subhorizontal discontinuity interfaces and shallow basins. J. Geophys. Res. 116, B02413. http://dx.doi.org/10.1029/2010JB008023. Camacho, A.G., Montesinos, F.G., Vieira, R., 2000. Gravity inversion by means of growing bodies. Geophysics 65, 95–101. Camacho, A., Montesinos, F., Vieira, R., Arnoso, J., 2001. Modelling of crustal anomalies of Lanzarote (Canary Islands) in light of gravity data. Geophys. J. Int. 147, 403–414. Canales, J., Dañobeitia, J., 1998. The Canary Islands swell: a coherence analysis of bathymetry and gravity. Geophys. J. Int. 132, 479–488. Christensen, N.I., Mooney, W.D., 1995. Seismic velocity structure and composition of the continental crust: a global view. J. Geophys. Res. 100, 9761–9788. Connolly, J., Kerrick, D., 2002. Metamorphic controls on seismic velocity of subducted oceanic crust at 100–250 km depth. Earth Planet. Sci. Lett. 204, 61–74.

J. Fullea et al. / Earth and Planetary Science Letters 409 (2015) 71–88

Connolly, J.A.D., 2005. Computation of phase equilibria by linear programming: a tool for geodynamic modeling and its application to subduction zone decarbonation. Earth Planet. Sci. Lett. 236, 524–541. Contrucci, I., Klingelhoefer, F., Perrot, J., Bartolomé, R., Gutscher, M.-A., Sahabi, M., Malod, J., Rehault, J.-P., 2004. The crustal structure of the NW Moroccan continental margin from wide-angle and reflection seismic data. Geophys. J. Int. 159, 117–128. Dañobeitia, J., Canales, J., 2000. Magmatic underplating in the Canary Archipelago. J. Volcanol. Geotherm. Res. 103, 27–41. Dañobeitia, J., Canales, J., Dehghani, G., 1994. An estimation of the elastic thickness of the lithosphere in the Canary Archipelago using admittance function. Geophys. Res. Lett. 21, 2649–2652. Dick, H., 1989. Abyssal peridotites, very slow spreading ridges and ocean ridge magmatism. Geol. Soc. (Lond.) Spec. Publ. 42, 71–105. Filmer, P., McNutt, M., 1989. Geoid anomalies over the Canary Islands group. Mar. Geophys. Res. 11, 77–87. Freundt, A., Schmincke, H.U., 1995. Petrogenesis of rhyolite–trachyte–basalt composite ignimbrite P1, Gran Canada, Canary Islands. J. Geophys. Res. 100, 455–474. Fullea, J., Afonso, J.C., Connolly, J.A.D., Fernàndez, M., García-Castellanos, D., Zeyen, H., 2009. LitMod3D: an interactive 3-D software to model the thermal, compositional, density, seismological, and rheological structure of the lithosphere and sublithospheric upper mantle. Geochem. Geophys. Geosyst. 10, Q08019. http://dx.doi.org/10.1029/2009GC002391. Fullea, J., Fernàndez, M., Zeyen, H., 2006. Lithospheric structure in the AtlanticMediterranean transition zone (southern Spain, northern Morocco): a simple approach from regional elevation and geoid data. C. R. Geosci. 338, 140–151. Fullea, J., Fernàndez, M., Zeyen, H., 2008. FA2BOUG – a FORTRAN 90 code to compute Bouguer gravity anomalies from gridded free air anomalies: application to the Atlantic-Mediterranean transition zone. Comput. Geosci. 34, 1665–1681. Fullea, J., Fernàndez, M., Zeyen, H., Vergés, J., 2007. A rapid method to map the crustal and lithospheric thickness using elevation, geoid anomaly and thermal analysis. Application to the Gibraltar Arc System and adjacent zones. Tectonophysics 430, 97–117. Fullea, J., Lebedev, S., Agius, M., Jones, A., Afonso, J.C., 2012. Lithospheric structure in the Baikal–central Mongolia region from integrated geophysical–petrological inversion of surface-wave data and topographic elevation. Geochem. Geophys. Geosyst. 13, Q0AK09. http://dx.doi.org/10.1029/2012GC004138. Fullea, J., Rodríguez-González, J., Charco, M., Martinec, Z., Negredo, A., Villaseñor, A., 2015. Perturbing effects of sub-lithospheric mass anomalies in GOCE gravity gradient and other gravity data modelling: application to the AtlanticMediterranean transition zone. Int. J. Appl. Earth Obs. Geoinf. 35, 54–69. http://dx.doi.org/10.1016/j.jag.2014.02.003. García-Castellanos, D., 2002. Interplay between lithospheric flexure and river transport in foreland basins. Basin Res. 14, 89–104. García-Yeguas, A., Koulakov, I., Ibáñez, J.M., Rietbrock, A., 2012. High resolution 3D P wave velocity structure beneath Tenerife Island (Canary Islands, Spain) based on tomographic inversion of active-source data. J. Geophys. Res. 117, B09309. http://dx.doi.org/10.1029/2011JB008970. Geldmacher, J., Hoernle, K., Bogaard, P.v.d., Duggen, S., Werner, R., 2005. New 40 Ar/39 Ar age and geochemical data from seamounts in the Canary and Madeira volcanic province: support for the mantle plume hypothesis. Earth Planet. Sci. Lett. 237, 85–101. González, P.J., Samsonov, S.V., Pepe, S., Tiampo, K.F., Tizzani, P., Casu, F., Fernández, J., Camacho, A.G., Sansosti, E., 2013. Magma storage and migration associated with the 2011–2012 El Hierro eruption: implications for crustal magmatic systems at oceanic Island volcanoes. J. Geophys. Res., Solid Earth 118, 4361–4377. http://dx.doi.org/10.1002/jgrb.50289. Gottsmann, J., Camacho, A., Martí, J., Wooller, L., Fernández, J., Garcia, A., Rymer, H., 2008. Shallow structure beneath the Central Volcanic Complex of Tenerife from new gravity data: implications for its evolution and recent reactivation. Phys. Earth Planet. Inter. 168, 212–230. Griffin, W.L., O’Reilly, S.Y., Ryan, C.G., 1999. The composition and origin of subcontinental lithospheric mantle. In: Fei, Y., Berkta, C.M., Mysen, B.O. (Eds.), Mantle Petrology: Field Observations and High-Pressure Experimentation: A Tribute to Francis R. (Joe) Boyd. In: Spec. Publ. – Geochem. Soc., vol. 6, pp. 13–45. Herzberg, C., 2011. Basalts as temperature probes of Earth’s mantle. Geology 39, 1179–1180. http://dx.doi.org/10.1130/focus122011.1. Herzberg, C., Asimow, P.D., 2008. Petrology of some oceanic island basalts: PRIMELT2.XLS software for primary magma calculation. Geochem. Geophys. Geosyst. 9, Q09001. http://dx.doi.org/10.1029/2008GC002057. Herzberg, C., Asimow, P.D., Arndt, N., Niu, Y., Lesher, C., Fitton, J., Cheadle, M., Saunders, A., 2007. Temperatures in ambient mantle and plumes: constraints from basalts, picrites, and komatiites. Geochem. Geophys. Geosyst. 8, Q02006. http://dx.doi.org/10.1029/2006GC001390. Hildenbrand, T., Kucks, R., Hamouda, M., Bellot, A., 1988. Bouguer gravity map and related filtered anomaly maps of Morocco. U.S. Geological Survey Open-File Report 88-517, 15 pp. Hoernle, K., Schmincke, H.-U., 1993. The role of partial melting in the 15-Ma geochemical evolution of Gran Canaria: a blob model for the Canary hotspot. J. Petrol. 34, 599–626.

87

Hoernle, K., Zhang, Y.-S., Graham, D., 1995. Seismic and geochemical evidence for large-scale mantle upwelling beneath the eastern Atlantic and western and central Europe. Nature 374, 34–39. Holland, T.J.B., Powell, R., 1998. An internally-consistent thermodynamic dataset for phases of petrological interest. J. Metamorph. Geol. 16, 309–344. Karato, S.i., 1993. Importance of anelasticity in the interpretation of seismic tomography. Geophys. Res. Lett. 20, 1623–1626. Klingelhoefer, F., Labails, C., Cosquer, E., Rouzo, S., Geli, L., Aslanian, D., Olivet, J.-L., Sahabi, M., Nouze, H., Unternehr, P., 2009. Crustal structure of the SW-Moroccan margin from wide-angle and reflection seismic data (the DAKHLA experiment) Part A: wide-angle seismic models. Tectonophysics 468, 63–82. Krastel, S., Schmincke, H.-U., 2002. Crustal structure of northern Gran Canaria, Canary Islands, deduced from active seismic tomography. J. Volcanol. Geotherm. Res. 115, 153–177. Le Bas, M., Rex, D., Stillman, C., 1986. The early magmatic chronology of Fuerteventura, Canary Islands. Geol. Mag. 123, 287–298. Legendre, C., Meier, T., Lebedev, S., Friederich, W., Viereck-Götte, L., 2012. A shear wave velocity model of the European upper mantle from automated inversion of seismic shear and surface waveforms. Geophys. J. Int. 191, 282–304. http://dx.doi.org/10.1111/j.1365-246X.2012.05613.x. Lodge, A., Nippress, S.E.J., Rietbrock, A., García-Yeguas, A., Ibáñez, J.M., 2012. Evidence for magmatic underplating and partial melt beneath the Canary Islands derived using teleseismic receiver functions. Phys. Earth Planet. Inter. 212, 44–54. Martinec, Z., Fullea, J., 2015. A refined model of sedimentary rock cover in the southeastern part of the Congo basin from GOCE gravity and vertical gravity gradient observations. Int. J. Appl. Earth Obs. Geoinf. 35, 70–87. http://dx.doi.org/10.1016/j.jag.2014.03.001. Martinez-Arevalo, C., Mancilla, F.D.L., Helffrich, G., Garcia, A., 2013. Seismic evidence of a regional sublithospheric low velocity layer beneath the Canary Islands. Tectonophysics 608, 586–599. McDonough, W.F., Sun, S.-S., 1995. The composition of the Earth. Chem. Geol. 120, 223–253. McKenzie, D., Bickle, M.J., 1988. The volume and composition of melt generated by extension of the lithosphere. J. Petrol. 29, 625–679. Minster, J.B., Anderson, D.L., 1981. A model of dislocation-controlled rheology for the mantle. Philos. Trans. R. Soc. A 299, 319–356. Montesinos, F., Arnoso, J., Benavent, M., Vieira, R., 2006. The crustal structure of El Hierro (Canary Islands) from 3-D gravity inversion. J. Volcanol. Geotherm. Res. 150, 283–299. Montesinos, F.G., 1999. Inversión gravimétrica 3D por técnicas de evolución. Aplicación a la isla de Fuerteventura. PhD thesis. Universidad Complutense de Madrid, Madrid. Neumann, Else-Ragnhild, 1991. Ultramafic and mafic xenoliths from Hierro, Canary Islands: evidence for melt infiltration in the upper mantle. Contrib. Mineral. Petrol. 106 (2), 236–252. Neumann, E.-R., Griffin, W.L., Pearson, N.J., O’Reilly, S.Y., 2004. The evolution of the upper mantle beneath the Canary Islands: information from trace elements and Sr isotope ratios in minerals in mantle xenoliths. J. Petrol. 45, 2573–2612. Neumann, E.-R., Wulff-Pedersen, E., Johnsen, K., Andersen, T., Krogh, E., 1995. Petrogenesis of spinel harzburgite and dunite suite xenoliths from Lanzarote, eastern Canary Islands: implications for the upper mantle. Lithos 35, 83–107. Neumann, E.-R., Wulff-Pedersen, E., Pearson, N., Spencer, E., 2002. Mantle xenoliths from Tenerife (Canary Islands): evidence for reactions between mantle peridotites and silicic carbonatite melts inducing Ca metasomatism. J. Petrol. 43, 825–857. Pavlis, N.K., Holmes, S.A., Kenyon, S.C., Factor, J.K., 2008. An earth gravitational model to degree 2160: EGM2008, EGU General Assembly 2008, 13–18. Pous, J., Heise, W., Schnegg, P.-A., Muñoz, G., Martí, J., Soriano, C., 2002. Magnetotelluric study of the Las Canadas caldera (Tenerife, Canary Islands): structural and hydrogeological implications. Earth Planet. Sci. Lett. 204, 249–263. Putirka, K.D., 2005. Mantle potential temperatures at Hawaii, Iceland, and the mid-ocean ridge system, as inferred from olivine phenocrysts: evidence for thermally driven mantle plumes. Geochem. Geophys. Geosyst. 6, Q05L08. http://dx.doi.org/10.1029/2005GC000915. Putirka, K., 2008. Excess temperatures at ocean Islands: implications for mantle layering and convection. Geology 36 (4), 283–286. Putirka, K.D., Perfit, M., Ryerson, F., Jackson, M.G., 2007. Ambient and excess mantle temperatures, olivine thermometry, and active vs. passive upwelling. Chem. Geol. 241, 177–206. Raffone, N., Chazot, G., Pin, C., Vannucci, R., Zanetti, A., 2009. Metasomatism in the lithospheric mantle beneath Middle Atlas (Morocco) and the origin of Fe- and Mg-rich wehrlites. J. Petrol. 50, 197–249. Ranero, C., Torne, M., Banda, E., 1995. Gravity and multichannel seismic reflection constraints on the lithospheric structure of the Canary Swell. Mar. Geophys. Res. 17, 519–534. Roest, W.R., Dañobeitia, J.J., Verhoef, J., Collette, B.J., 1992. Magnetic anomalies in the Canary Basin and the mesozoic evolution of the Central North Atlantic. Mar. Geophys. Res. 14, 1–24.

88

J. Fullea et al. / Earth and Planetary Science Letters 409 (2015) 71–88

Shibata, T., Thompson, G., 1986. Peridotites from the Mid-Atlantic Ridge at 43◦ N and their petrogenetic relation to abyssal tholeiites. Contrib. Mineral. Petrol. 93, 144–159. Simon, N.S., Neumann, E.-R., Bonadiman, C., Coltorti, M., Delpech, G., Grégoire, M., Widom, E., 2008. Ultra-refractory domains in the oceanic mantle lithosphere sampled as mantle xenoliths at ocean islands. J. Petrol. 49, 1223–1251. Smith, W.H.F., Sandwell, D.T., 1994. Bathymetry prediction from dense satellite altimetry and sparse shipborne bathymetry. J. Geophys. Res. 99, 21803–21824. Smith, W.H.F., Sandwell, D.T., 1997. Global seafloor topography from satellite altimetry and ship depth soundings. Science 277, 195–196. Spieker, K., Wölbern, I., Thomas, C., Harnafi, M., El Moudnib, L., 2014. Crustal and upper-mantle structure beneath the western Atlas Mountains in SW Morocco derived from receiver functions. Geophys. J. Int. 198 (3), 1474–1485. Verhoef, J., Collette, B.J., Dañobeitia, J.J., Roeser, H.A., Roest, W.R., 1991. Magnetic anomalies off west-Africa (20–38◦ N). Mar. Geophys. Res. 13, 81–103. Vieira, R., Toro, C., Araña, V., 1986. Microgravimetric survey in the caldera of Teide, Tenerife, Canary Islands. Tectonophysics 130, 249–257.

Watts, A., 1994. Crustal structure, gravity anomalies and flexure of the lithosphere in the vicinity of the Canary Islands. Geophys. J. Int. 119, 648–666. Watts, A.B., 2001. Isostasy and Flexure of the Lithosphere. Cambridge University Press. Watts, A., Peirce, C., Collier, J., Dalwood, R., Canales, J., Henstock, T., 1997. A seismic study of lithospheric flexure in the vicinity of Tenerife, Canary Islands. Earth Planet. Sci. Lett. 146, 431–447. White, W.M., 2010. Oceanic island basalts and mantle plumes: the geochemical perspective. Annu. Rev. Earth Planet. Sci. 38, 133–160. Wittig, N., Pearson, D.G., Baker, J.A., Duggen, S., Hoernle, K., 2010. A major element, PGE and Re–Os isotope study of Middle Atlas (Morocco) peridotite xenoliths: evidence for coupled introduction of metasomatic sulphides and clinopyroxene. Lithos 115, 15–26. Wulff-Pedersen, E., Neumann, E.-R., Jensen, B., 1996. The upper mantle under La Palma, Canary Islands: formation of Si–K–Na-rich melt and its importance as a metasomatic agent. Contrib. Mineral. Petrol. 125, 113–139. Ye, S., Canales, J., Rihm, R., Danobeitia, J., Gallart, J., 1999. A crustal transect through the northern and northeastern part of the volcanic edifice of Gran Canaria, Canary Islands. J. Geodyn. 28, 3–26.