Anelastic response of the Earth's crust underneath the Canary Islands revealed from ocean tide loading observations

SUMMARY

We report on the analysis of M₂ ocean tide loading (OTL) kinematic GPS vertical displacement and tidal gravity measurements using 26 GPS and four gravimetric sites across the Canary Islands archipelago. In this region, the standard deviation among recent ocean tide models is lower than 0.4 cm in amplitude and 0.3° in phase, which are suitably accurate for displacement modelling. However, for gravity we need to model regional ocean tides to achieve enough accuracy in the loading calculations. Particularly, this study improves the predicted OTL gravity variations when global ocean models are replaced with the regional model CIAM2 which assimilates local tide gauge data. These small ocean tide model errors allow us to use the differences between observed and predicted OTL values to study the elastic and anelastic properties of the solid Earth around the Canary Islands. In the prediction of OTL, we first used the recent elastic STW105 and S362ANI seismic models, obtaining average observed minus predicted residuals of 1.2–1.3 mm for vertical displacement and 3 nm s⁻² for gravity. After the STW105 and S362ANI models were adjusted for anelasticity, by considering a constant quality factor Q at periods ranging from 1 s to 12.42 hr, the average misfit between observations and predicted OTL values reduced to 0.7–0.8 mm for vertical displacement and to 1 nm s⁻² for gravity. However, the average vertical displacement misfit is made up from site misfits less than 0.5 mm in western islands but for the easternmost islands of Lanzarote and Fuerteventura, they still reach up to nearly 2 mm at some sites, which still exceeds the uncertainty in the GPS observations. It is hypothesized that mantle upwelling underneath the Canary Islands, creating spatial variations in the elastic properties, causes the large residuals observed in the eastern islands. We reduced the shear modulus by up to 35 per cent in the upper mantle layer of 24.4–220 km depth. This produced residual observed minus model differences of about 0.7 mm for the sites on Lanzarote and Fuerteventura, comparable to the results obtained for the GPS sites across the rest of the archipelago, whose residuals in turn were also slightly reduced through the V_S velocity and shear modulus reductions (by 0.2 mm on average).

1 INTRODUCTION

The weight of the ocean tides produces a periodic load, called ocean tide loading (OTL), that deforms the Earth's surface and causes gravity changes. OTL can cause vertical displacements of several centimetres on the Earth's surface that can be measured with Global Positioning System (GPS) receivers (e.g. Allinson et al. 2004; Penna et al. 2015; Martens et al. 2016). Similarly, OTL produces gravity variations of a few hundred nm s⁻² in amplitude, which can be measured with carefully calibrated relative or absolute gravimeters (e.g. Farrell 1972; Baker & Bos 2003). If models of the ocean tides are accurate, the magnitude of these OTL deformations can be used to provide information about the interior structure of the solid Earth (Farrell 1973).

Until about the last decade, errors in ocean tide models were the limiting factor in the prediction of OTL values, although the elastic properties of the Earth used to compute OTL were based on global seismic models such as the Preliminary Reference Earth Model (PREM, Dziewonski & Anderson 1981). However, ocean tide models have now become more accurate and the limitations of using PREM to compute OTL have become apparent (e.g. Ito & Simons 2011; Bos et al. 2015; Martens & Simons 2020; Wang et al. 2020). Seismic earth models such as PREM are valid at the reference period of 1 s, whereas at tidal periods the elastic properties are slightly different due to anelastic dispersion (Liu et al. 1976; Benjamin et al. 2006), which principally happens in the asthenosphere. If the OTL signal is large enough, say larger than 20 mm in amplitude for vertical displacements, and if the stations are within a 250 km distance from the coast, then this effect is large enough to be observed in tidal GPS and gravity observations.

The tidal periods fall between the seismic periods (around 1 s) and the very much slower processes such as mantle convection and flexure of the crust, processes that take place over tens of millions of years. Although at seismic periods the upper mantle of the solid Earth behaves as a solid, at longer timescales the material behaves more viscously or even as a fluid (e.g. in the asthenosphere). However, how the elastic properties of the upper mantle change from periods of 1 s to 10 Myr is still not well understood (Lau & Holtzman 2019). Liu et al. (1976) introduced a power-law model for the quality factor Q, which is inversely proportional to the amount of dissipation during a cycle, with an exponent α that depends on the physical process causing the dissipation. This Q^α model describes the attenuation of the body and surface waves and the free oscillation of the Earth. Zschau (1978) was one of the first to use a Q^α model to modify the elastic constants of seismic earth models for OTL computations. A general model for frequency dependence of Q across the seismic band is in the form (Benjamin et al. 2006)

where ω is the frequency, Q₀ is the value of Q at some reference seismic frequency ω₀. The shear modulus can be written as |$\mu \ = {\mu }_0{\rm{\ }} + {\rm d}\mu $|⁠, where |${\mu }_0$| is the shear modulus at the reference period of 1 s (at |${\omega }_0$|⁠) and |${\rm d}\mu $| is the change in its value due to anelastic dispersion at frequency |$\omega $|⁠. It can be computed as follows (Benjamin et al. 2006)

Thus, we introduce complex-valued shear moduli, as in Zschau (1978). Since the quality factor for the bulk modulus is very high, we can assume that it is not affected by anelastic dispersion, which provides us with the following relation for the change in the first Lamé parameter λ

With these modified elastic constants, one can compute new complex-valued load Love numbers (Farrell 1972; Bos & Scherneck 2013) and from this set of load Love numbers, new Green's functions may be formed. Following Bos et al. (2015) only the real part of the Green's functions has been used.

Despite the theoretical framework for modelling the change in elastic properties of the solid Earth over a wide range of frequencies being well advanced, more high accuracy observations are needed to constrain the parameters of these models. Previous works on this topic are those of Ito & Simons (2011) and Yuan & Chao (2012) who observed OTL for the western United States using GPS. Bos et al. (2015) reported the importance of anelasticity in M₂ OTL computations for GPS vertical displacements in western Europe, which is a passive margin zone. More recently, Wang et al. (2020) investigated OTL around the East China Sea region (which overlies a subduction zone). All these studies found indications that the asthenosphere is less stiff at tidal periods.

In this paper, the study area is an intraplate volcanic archipelago lying on oceanic lithosphere, which has different tectonic characteristics than the aforementioned works. Also, the elastic and anelastic properties of the upper layers (depth less than 200 km) of the solid Earth are adjusted to represent better the average conditions underneath the Canary Islands. The accuracy of these adjustments will be validated with tidal kinematic GPS vertical displacements and new tidal gravity data. The calculation involved in this validation uses recent global ocean tide models, each combined with a regional one that is based on local tide gauge and satellite altimetry data assimilated into a hydrodynamic model. The part of the global model covering the regional model's area was removed and replaced with the regional model, to create a hybrid (hereafter we refer to this as ‘global model replaced with regional one’). Differences between the observed and predicted/computed tidal deformation provide information about the elastic properties of the solid Earth, since such differences are not dominated by errors in the ocean tide model. In addition, we use earth models adjusted for anelastic dispersion from the reference period of 1 s to the period of 12.42 hr (harmonic M₂), which mostly affects the asthenosphere. We study the discrepancies between GPS vertical displacements and gravity measurements and the respective OTL predictions to investigate some anelastic effects in the upper mantle structure below this ocean island volcanic archipelago. We consider only the harmonic M₂ due to its relatively large amplitude over the Canary Islands region.

2 GEODYNAMIC FRAMEWORK

The Canary Islands archipelago is situated in the Atlantic Ocean close to the Moroccan passive continental margin, to the east of the Central Atlantic Ocean (Fig. 1). The islands Lanzarote and Fuerteventura are along a north–northeast-trending submarine ridge (the East Canary Ridge) that extends for about 250 km parallel to the African coast and may be aligned along the contact between attenuated continental crust on the east and oceanic crust on the west. The archipelago is located on the African Plate and lies on Jurassic (150–170 Ma) oceanic lithosphere. Seismic refraction experiments and geomagnetic data pointed out that the transition between oceanic crust to continental crust is located around Lanzarote, Fuerteventura and the northwest coast of Africa. The islands form independent volcanic edifices separated by narrow channels, except for Lanzarote and Fuerteventura that share the aforementioned submarine ridge (e.g. Schmincke 1982; Roest et al. 1992; Ancochea et al. 1996, 2004; Watts et al. 1997; Carracedo 1999; Acosta et al. 2005 ).

The Canary Islands were created from volcanic and tectonic activity that started around 70–80 Ma. The oldest subaerial volcanism located in the eastern islands (Fuerteventura and Lanzarote from about 20 Ma) and the youngest one in the westernmost islands (El Hierro and La Palma less than 2 Ma) conforms with an east–west age decrease trend. However, the lack of regularity in this age-decrease, the existence of several Ma gaps in the volcanic activity, the irregular distribution of islands and seamounts in the Canary province, and the lack of a broad topographic swell as observed in Hawaii and Cape Verde archipelagos, cannot be explained by the conventional hotspot model, that is, the islands were formed as the African plate moved over a hotspot (e.g. see Guillou et al. 1996; Filmer & McNutt 1989; Ye et al. 1999; Geldmacher et al. 2005; Lodge et al. 2012; van der Bogaard 2013). The unifying model proposed by Anguita & Hernán (2000), which is the most accepted model to explain the genesis of the archipelago, supports the existence of a residual plume under north Africa, the Canary Islands and western and central Europe and a fracture system that allows magma to ascend. Therefore, the volcanic edifices that comprise the archipelago resulted from the mantle plume upwelling whose dynamics produced lithospheric melting periods of active volcanism from the Jurassic to present (van der Bogaard et al. 2013). Duggen et al. (2009) suggested a mantle plume located to the west of the archipelago, which facilitates the flow of the mantle material to the east, and crossing the northwest of Africa. The mantle material flows through the base of the oceanic lithosphere and into a subcontinental lithospheric corridor that was already inferred from geophysical experiments (Fullea Urchulutegui et al. 2006; Missenard et al. 2006; Teixell et al. 2005). These authors concluded that there is no conflict between the mantle flow model and the involvement of lithospheric processes in magma generation and ascent in the Canary Islands. Thus, the kinematics of the Eurasian and African plates, the Mid-Atlantic Ridge and other regional structures control the plate-scale stress field, and together with the regional fault systems influence the tectonics in the archipelago (Geyer et al. 2016). Furthermore, recent regional GPS studies (Arnoso et al. 2020; Barbero et al. 2021) confirm areas of high strain rates (mostly NW–SE and NE–SW trending) associated with extensional and compressional tectonic regimes controlling magma accumulation and episodes of volcanism in that region.

In this geological context, Lanzarote and Fuerteventura are located in the emergent part of the East Canary Ridge, to the eastern margin of the archipelago, and close to continental Africa (e.g. Marinoni & Pasquaré 1994). This ridge consists of several uplifted blocks of oceanic basement, with a sedimentary cover more than 10 km thick (Banda et al. 1981). Considering the aforementioned unifying approach of Anguita & Hernan (2000), Duggen et al. (2009) explain that thinning of the lithosphere to the east of the archipelago enhances the mantle flow along its base to areas of thinner lithosphere, since Neumann et al. (1995), based on petrological data, explained that the lithosphere beneath the easternmost Canary island is thinner than beneath the more western islands. Dañobeitia & Canales (2000) supported a 7–12-km-thick layer, interpreted as magmatic underplating at the eastern islands, and suggested that the interaction of the mantle plume with the oceanic crust could be responsible for the crustal thickening observed beneath the islands, which can be explained by their proximity to the African margin. More recently, Lodge et al. (2012) determined for Lanzarote a velocity model with discontinuities at approximately 4, 10 and 18 ± 1.5 km depth (the Moho boundary), suggesting either a magmatic layer of underplating in the oceanic crust up to 8 km beneath this island, or that Lanzarote is on a region of thinned continental crust. There are geophysical studies which consider the complex crustal structure of the lithosphere in the Canary Islands (Banda et al. 1981; Watts et al. 1997; Ye et al. 1999; Dañobeitia & Canales 2000; Contrucci et al. 2004), reporting different discontinuities in depth, in particular, to the east of the archipelago in the vicinity of the Moroccan continental margin. Thus, the Moroccan passive margin to the north of the Canary Islands is characterized by a progressive transition from old oceanic crust to a 34–35-km-thick continental crust, whereas to the south of the Canary Islands the continental margin is 27-km-thick crust, and the Moho boundary is at about 27 km depth (Klingelhoefer et al. 2009). To the west of the Moroccan Atlas and close to the Atlantic margin, Spieker et al. (2014) estimated the Moho depth as 24.65 to 28.05 ± 6.7 km. The geophysical-petrological approach of Fullea et al. (2015) points to a thin lithosphere about 100 km thick overlying an anomalous asthenosphere 100 °C hotter than the normal sublithospheric mantle beneath the archipelago. Finally, the recent study of Negredo et al. (2022) modelled the ascent of an upper mantle plume resulting in an asymmetric mantle flow pattern and lateral migration due to the interaction with the heterogeneous lithosphere.

3 ANALYSIS OF GPS AND GRAVITY TIME-SERIES

The number of accessible continuous GPS public sites across the Canary Islands has increased recently due to the 2004–2005 seismic crisis detected on Tenerife and the 2011–2012 submarine eruption in El Hierro, which revealed the need for improved ground deformation studies in the archipelago. We collated GPS data from 26 sites logging over the time window 2009.0–2016.8 (Fig. 1 and Table 1) from several data stores (see Data Availability section). Data durations per GPS site range from about 1–7.8 yr (the median time-series length is nearly 6 yr), as listed in Table 1.

The GPS data from all sites were analysed in kinematic precise point positioning (PPP) mode using the NASA GIPSY v6.4 software, applying the optimized processing settings for OTL displacement measurement determined by Penna et al. (2015) using data from western Europe. Their tests, based on the insertion and recovery of synthetic harmonic displacement signals, suggested a vertical OTL displacement measurement accuracy of better than 0.5 mm. Processing options included Jet Propulsion Laboratory (JPL) repro2 GPS fixed fiducial orbits and 30 second clocks, European Centre for Medium-Range Weather Forecasts (ECMWF) model zenith hydrostatic and wet delays (Boehm et al. 2006), with the zenith wet delay estimated via a random walk process [0.1 mm (√s)⁻¹ process noise] every 5 min together with the site coordinates [3.2 mm (√s)⁻¹ process noise]. A carrier phase observational standard deviation of 10 mm was used, and the VMF1 gridded tropospheric mapping function was applied. Earth body and pole tides were modelled following the International Earth Rotation Service Conventions (Petit & Luzum 2010) and OTL predicted corrections were applied using the Gutenberg–Bullen earth model and the FES2004 ocean tide model (Lyard et al. 2006), computed using the SPOTL v3.2.4 program (Agnew 1997). Wang et al. (2020) used the same GIPSY GPS processing settings as Penna et al. (2015) to estimate M₂ OTL displacement over the East China Sea, having first checked that a synthetic harmonic signal could also be recovered to better than 0.5 mm, with Abbaszadeh et al. (2020) also finding similar optimal process noise settings for 41 GPS sites globally using the PANDA software (and validated against GIPSY).

For each GPS site, the 5-min PPP coordinate estimates were averaged to one value every 30 min and then time-series formed by concatenation. Then, the tidal parameters for the GPS vertical component were estimated via harmonic analysis with the ETERNA software (Wenzel 1996), with an average amplitude standard deviation of ±0.3 mm. At the M₂ frequency, the GPS-estimated vertical displacement residuals were added to the predicted OTL values to obtain the total GPS-observed OTL displacement per site, as per Bos et al. (2015).

The gravity observations used in this study were obtained with the Micro-g LaCoste & Romberg spring gravimeter gPhone-054 model on El Hierro (AULA site) and Tenerife (IZA site), and the Graviton EG-1194 model on Gran Canaria (CVAR site) and Lanzarote (TIM site; Fig. 1). The performance of these meters in terms of resolution, accuracy and long-term stability was studied when compared to the superconducting gravimeters SG-C026 and OSG-064, at sites J9 in Strasbourg (France) and CDT-Yebes (Spain), respectively (Riccardi et al. 2011a; Arnoso et al. 2014). Those comparisons provide a definition of the spring gravimeter's transfer function in the tidal band, in terms of phase and amplitude. Also, the scale factors of both instruments, which showed variations of about 0.1 per cent, were determined by colocated measurements with absolute gravimeters at the superconducting gravity observatories. In addition, Rosat et al. (2015) studied the respective operation of those instruments in terms of noise levels through the computation of the residual power spectral densities, and found similar noise magnitudes for both spring gravimeters at the subseismic frequency band. The tidal data were all analysed with the VAV software (Venedikov et al. 2003, 2005), which is based on the tidal harmonic analysis method. VAV has the capacity to find anomalous data in the processing through comparison of the gravity residuals with a fixed threshold level of significance for a selected probability confidence interval. This is done in a recursive procedure where the data having residuals equal to or higher than the threshold are eliminated from the processing and, accordingly, the precision of the estimation of tidal parameters increases (see e.g. Arnoso et al. 2011).

4 EARTH CRUST MODELS, ELASTIC GREEN'S FUNCTIONS AND OCEAN TIDE LOADING CALCULATION

4.1. Earth models and elastic Green's functions

For our earth models, we have used two recent global models: STW105 and S362ANI (Kustowski et al. 2008); since the latter varies spatially, we have used the average depth profile beneath the Canary Islands. Fig. 2(a) shows depth profiles of density (ρ), P-wave (V_P) and S-wave (V_S) velocities for these two models, as well as for PREM. STW105 and S362ANI only differ for V_S. Since we found that the loads computed using PREM show systematically poorer fits we discuss them only in the Supporting Information (Section S1). We also examined models that use the crustal structure for Fuerteventura, Gran Canaria and Lanzarote (Fig. S1). Because this changes the Green's functions only for distances less than 10 km, the differences in OTL calculation between these models and STW105 are only 0.2–0.4 mm and 0.7–1.0 nm s⁻² for M₂ vertical displacement and gravity, respectively (Fig. S2 and Table S1). This is small enough that we did not use these local crustal models.

As mentioned above, the elastic properties of the Earth's crust could be different at lower frequencies than the seismic reference period of 1 s due to anelastic dispersion effects in the asthenosphere. They are derived for a reference period of 1 s and adjusted to be valid at harmonic M₂ (12.42 hr), assuming that anelastic dispersion takes place in the asthenosphere with a constant quality factor Q (Bos et al. 2015). We only modified the shear modulus because Q in the asthenosphere for the bulk modulus is substantially larger and therefore less influenced by anelastic dispersion. The depth profiles for the quality factor Q, as well as for the bulk and shear modulus, are shown in Figs 2(b)–(d), respectively. These plots include the variation applied when lowering the quality factor Q to 75, 50 and 35, assuming a Q^α model with α = 0.1. For our calculations we used in each layer the Q values tabulated by each earth model.

Fig. 3(a) displays the Green's functions for the models STW105 and S362ANI, both modified for anelastic dispersion effects, and the elastic version of STW105 for comparison. Fig. 3(b) shows the ratios of the Green's functions after dividing models STW105 and S362ANI, corrected for anelastic dispersion, by the elastic STW105. It can be observed that the change mainly affects the Green's function for distances of 10–100 km from the calculation point. Also, Fig. 3 includes plots of the Green's functions for the modified STW105 model after lowering the constant quality factor Q. In Section 5.2 we present a sensitivity study regarding this variation of Q for the asthenosphere.

4.2 Ocean tide models

We used the global ocean tide models FES2014b (Lyard et al. 2021), GOT4.10c (Ray 2013) and TPXO9v5a (Egbert et al. 2002). These models (and their respective previous versions) have been tested by several authors, who found standard deviations among the models for M₂ tidal elevations of about 0.3 cm for the deep oceans and up to 3.4 cm over shallow water areas (e.g. Zaron & Elipot 2020; Stammer et al. 2014). Additionally, we studied these global models in the Canary Islands area bounded by coordinates 26.5°N to 30°N and 12.5°W to 19°W, and calculated the differences on a 0.0611° × 0.0611° grid, which yielded standard deviations of the differences lower than 0.4 cm in amplitude and 0.3° in phase. We also considered a regional ocean tide model specifically for the Canary Islands archipelago, called CIAM2 (Benavent 2011), which replaced the matching area of the three global models in order to more accurately retrieve the OTL. CIAM2 covers the aforementioned geographic area, and results from assimilating satellite altimetry and tide gauge data on several Canary Islands into a hydrodynamic model. This regional model yields a high accuracy ocean tide representation, a suitable coastline definition and water depth modelling. CIAM2 improved the previous CIAM1 model (Arnoso et al. 2006) by assimilating additional tide gauge data in the Canary Islands region. An extensive study of CIAM2 errors was carried out in Benavent (2011) by comparing CIAM2 with global and regional ocean tide models, and also through the Monte Carlo simulation approach (Dushaw et al. 1997). Basically, this simulation uses an iterative generation of inverse models dealing with synthetic data from hydrodynamic models (direct solutions) obtained, in turn, by random perturbation of the model parameters. The random errors, representing the noise and other non-tidal signals, are then added to the synthetic data. Once the process finishes, the difference between the direct and inverse solutions provides the estimation of the errors in the tidal heights. The simulation method gave errors not greater than 2 cm for the M₂ dominant wave along the CIAM2 domain. Section S3 of the Supporting Information includes a comparison (Tables S2 and S3 and Fig. S3) between the tidal amplitudes and phases, for harmonic M₂, observed by tide gauges and by satellite altimetry at the crossover points in the Canary Islands region with those amplitude and phases provided by ocean tide models (FES2014b, GOT4.10c, TPXO9v5a and CIAM2). Although the use of CIAM2 to replace the global ocean models locally does not influence the OTL vertical displacement calculated for any station, it improves the OTL gravity calculation in all stations used in this work.

4.3 Ocean tide loading calculation

Vertical displacements and gravity variations due to the OTL effect were predicted by convolving over the entire surface of the oceans (Ω) the tide distribution with the Green's functions (Fig. 3), G, calculated from a radially symmetric earth model, following the method of Farrell (1972). The OTL L, can be written in the space domain, through the integral

where ρ is the density of seawater, H is the complex-valued height of the ocean tide, r and r′ are the position vectors of the points of observation and of the loading, respectively.

The accuracy of the OTL values therefore depends on the accuracy of the ocean tide models and the Green's function that represents the elastic properties of the solid Earth. Furthermore, to compute the mass of the ocean tides that is pressing on the ocean floor, the seawater density is needed, which is assumed to be a constant with value 1030 kg m⁻³. The coastline also needs to be accurately modelled, especially for the gravity OTL values which are more sensitive to the very nearby ocean tides. All OTL computations were carried out using the software of Benavent (2011). Comparisons with the CARGA software (Bos & Baker 2005) produced similar results (for M₂ gravity, the difference is lower than 0.1 nm s⁻² in amplitude and 0.03° in phase) when using the same computing options, for example ocean tide model, Green's function, ocean-land mask for grid refinement.

For the OTL calculations, the oceanic grid on which the models are represented was recursively refined around the calculation point in order to more accurately fit the coastline. For a specific refinement area (in our case, the grid cells within 0.75° from the calculation point), each oceanic cell was subdivided into four smaller subcells until the distance to the calculation point was twenty times greater than the cell radius. High-resolution land–water mask files (cell size not smaller than 180 m) were used to determine the shorelines and, consequently, the land–water condition of the resulting refined cells. Thus, we built the land–water mask from the GEBCO data set (https://www.gebco.net) and digital terrain models from the National Geographic Institute of Spain (https://www.ign.es). Finally, we determined the tidal heights either by bilinear interpolation or by averaging of the surrounding cells from the ocean tide model.

Fig. 4 shows the OTL prediction calculated following the aforementioned procedure for harmonic M₂ across the Canary Islands using FES2014b replaced with CIAM2 and Green's functions based on the STW105 earth model with anelastic dispersion, both for vertical displacement and gravity. In accordance with M₂ ocean tide amplitudes in the archipelago, the amplitude of the OTL vertical displacement (Fig. 4a) is greater to the north of the archipelago, reaching the greatest value on the northeast of Lanzarote where it exceeds 2.9 cm, and reducing to about 2.3 cm on El Hierro in the southwest of the archipelago. Regarding gravity, the OTL amplitudes (Fig. 4b) are greater in the western islands (it can reach a value of about 130 nm s⁻² on La Palma and El Hierro) than in the central and eastern islands. The territory of the western islands is smaller in size than the other islands, however the altitude above sea level reaches up to about 2400 m (La Palma) or 1500 m (El Hierro, La Gomera) and the distances to the shoreline are short. More specifically, the respective ratios of elevation to distance for the four gravity stations are 0.33 (AULA), 0.20 (IZA), 0.07 (CVAR) and 0.06 (TIM). This fact influences the Newtonian component of the OTL calculation, which becomes similar in magnitude to the elastic component for those islands and, thus, increases the gravity OTL amplitudes. The OTL phases for both gravity and vertical displacement present a southwest–northeast trend, according to cophase maps of OTL models for M₂ around the Canary Islands (e.g. Arnoso et al. 2006), the highest (absolute) values being located in the eastern islands.

We also investigated the influence of ocean tide model errors considered in this study in the OTL computation (see Figs S4 to S6). Since the errors in the ocean tide models are unknown, we used the difference between the ocean tide models as a proxy. For harmonic M₂, we found that the differences in OTL vertical displacement between FES2014b and GOT4.10c or TPXO9v5a are smaller than 0.3 mm for the Canary Islands (Fig. S6 and Table S6). We checked the impact on the OTL calculation of each global model with and without CIAM2 replaced locally (Fig. S5 and Table S5), and we obtained mean values lower than −0.15 ± 0.33 mm for vertical displacement and 0.53 ± 2.46 nm s⁻² for gravity. From Figs S5 and S6 (and Tables S5 and S6) we can argue that, for OTL vertical displacement, none of the three global ocean models give better results than other, with or without using the regional model CIAM2. Hence for displacement, hereafter we only show results using OTL predicted with FES2014b replaced with CIAM2. For gravity, we found that the OTL differences between the three global models can reach up to 2.3 per cent of the amplitude. After replacing the global models with CIAM2 these differences reduce to about 0.3 per cent of the OTL amplitude (see statistics in Table S6). Furthermore, we will show in Section 5 that the variations in the amplitude of the gravity residuals at our four observation sites are within 0.2 nm s⁻² when using the three global models replaced with CIAM2.

5 RESULTS AND DISCUSSION

5.1 Tidal gravity residuals

Table 2 lists the tidal amplitudes and phases obtained for harmonic M₂ at the four observing sites, following the procedure described in Section 3. The lowest values of mean square deviations of the tidal amplitudes are obtained at AULA and IZA, both operated with gPhone gravimeters. The improved gPhone precisions over the Graviton model are attributed to the thermal system of the gPhone increasing the temperature stability, and its vacuum sealing making it almost insensitive to buoyancy changes due to atmospheric pressure fluctuations (Riccardi et al. 2011b). Both instruments are spring type gravimeters, whereby the instrumental drift affecting them influences the precise determination of the tidal parameters, if it is not properly modelled. To provide an idea of the noise around the tidal frequencies, we calculated the power spectra of the residuals for the four gravity sites and also for the GPS sites on Tenerife. Fig. S9 shows the tidal signal strongly reduced at all gravity stations and due to the large signal at M₂, there is some residual signal left at the semi-diurnal frequency. For GPS residuals time-series, the large residual signal at the daily period can be produced by the fact that we are analysing daily segments of GPS data. The daily signal also creates some false signals at the semidiurnal period. However, since this false signal seems to be random, its effect on the tidal analysis is limited (Penna et al. 2015).

The tidal gravity measurements were compared with the Dehant et al. (1999) body tide model. This model provides the analytical expressions used to compute numerical values of the tidal gravimetric factors for two earth models, one of them assuming an elastic Earth, and the other one including inelasticity and non-hydrostatic structure of the Earth. The difference between both models is of the order 0.12 per cent of the tidal amplitude, and less than 0.005° for the phase-delay (Dehant & Zschau 1989). The anelastic behaviour at tidal periods is still not well known throughout the Earth, and the residuals obtained from the comparison of the tidal gravity observations with the respective OTL predictions (calculated using Farrell's procedure described in Section 4.1) might enable the determination of the elastic or inelastic Earth's response in the Canary Islands, which can be used to study the effects of lateral inhomogeneities and/or mantle anelasticity in that region.

The phasor plots of Fig. 5 show the respective observed tidal gravity residuals for harmonic M₂, computed by subtracting the Dehant et al. (1999) body tide, considering both the elastic (DDWe) and inelastic (DDWi) models, from the observed values. Also plotted are the gravity OTL predictions using the three global ocean tide models replaced with the regional model CIAM2, and using Green's functions based on each of the S362ANI and STW105 earth models, all with anelastic dispersion effects. The misfit is about 1 nm s⁻² and models S362ANI and STW105 appear grouped and closer to the anelastic DDWi model, that is we find that the gravity residuals support an anelastic behaviour of the upper crust layers model at all four gravity sites along the archipelago.

The impact on the observed tidal gravity residuals and the gravity OTL of including the regional model CIAM2 in the OTL computation, in addition to the global FES2014b model, is shown in Fig. 6. It is evident that for all the gravity sites the accuracy in the computations is significantly higher when the CIAM2 model is also incorporated. Using CIAM2 and FES2014b, the amplitudes and phases of the OTL lead to a respective discrepancy of 0.0–0.5 per cent and −0.2° to 0.2° with respect to the observed tidal gravity residuals based on the inelastic DDWi non-hydrostatic earth model. However, using FES2014b alone without CIAM2, the respective discrepancies in amplitudes and phases are about −2.2 to 3.5 per cent and −0.4° to 0.6°. However, interchanging FES2014b for either GOT4.10c or TPXO9v5 and similarly replacing with CIAM2 does not substantially impact the result of the OTL calculation: arising maximum differences in the OTL gravity change for M₂ are in the range of −0.2 to 0.2 nm s⁻² across the archipelago (see Table S6 and Fig. S6).

5.2 GPS vertical displacements

Fig. S7(a) displays the vectors of the OTL response for GPS-observed and predicted vertical displacement for harmonic M₂ along the Canary archipelago, with the prediction calculated using the earth model STW105 with anelastic dispersion and the ocean tide model FES2014b replaced with the regional model CIAM2. The amplitude of the observed displacement ranges from 23.9 mm on El Hierro in the southwest of the archipelago to 29.3 mm on Lanzarote, and the mean value of the observed Greenwich phase lags is about −142.2°. For vertical displacement, OTL results vary by at most 0.2 mm using one or the other of the three global ocean tide models replaced with CIAM2, at our observation sites, so the choice of global model when replaced with CIAM2 does not matter. Fig. 7(a) displays the vector differences between the observed and predicted GPS vertical displacements, considering the Green's functions based on earth model STW105 with and without anelastic dispersion effects. It can be seen that the residuals are lower for all sites when the earth model includes the effect of the velocity dispersion in the asthenosphere at tidal periods, with the mean residual reducing from 1.3 to 0.8 mm (Table 3), and the greatest reduction of 0.5 mm at site TIAS (Lanzarote). It is worth noting that if the model FES2014b is not replaced with CIAM2 (see Figure S7b), using the earth model STW105 with anelastic dispersion, the result does not vary the interpretation of the vector differences, as the average value at the GPS sites with and without using CIAM2 differs by 0.1 mm (the maximum difference is about 0.4 mm). Similarly, for the earth model S362ANI with anelastic dispersion, we obtained the same values when using FES2014b only or FES2014b replaced with CIAM2. We could have used any of the global models FES2014b, TPXO9v5a or GOT4.10c without replacement with CIAM2 and it would not have impacted the displacement results substantially.

Next, we performed the comparison between observed GPS vertical displacements and predictions using different earth models (S362ANI and STW105, all including anelastic dispersion effects), all using FES2014b replaced with the regional ocean tide model CIAM2 as above. The respective residuals for each earth model are shown in Fig. 7(b) and listed numerically in Table 3. Regarding the earth models, it can be seen that the average misfit between observations and predicted OTL values is about 0.7–0.8 mm. The largest observed minus model residuals occur on Lanzarote and Fuerteventura, where the average magnitude of the residuals calculated for all GPS sites is about 1.5 mm, with residuals of 1.1–1.8 mm when using S362ANI and 1.3–1.9 mm when using STW105.

The remaining residuals (Fig. 7b) could represent noise in the observations or point to some effect of lateral heterogeneities in the Earth's structure, as was found by Martens & Simons (2020) in Alaska. The study of Latychev et al. (2009) proposes that the 3-D mantle structure can perturb the radial displacement and surface gravity within the semidiurnal band by ∼1.0 mm and ∼2 nm s⁻², respectively. We have obtained average GPS M₂ residuals of 0.7–0.8 mm when using earth models with a constant Q from 1 s to the period of M₂ in the asthenosphere, but the residuals at some sites on Lanzarote and Fuerteventura still reach up to 2 mm (Fig. 7b and Table 3). It is hypothesized that mantle upwelling underneath the Canary Islands, creating spatial variations in the elastic properties, causes the large residuals observed in the eastern islands. Studies from Teixell et al. (2003), Duggen et al. (2009) and, more recently, Miller et al. (2015), Civiero et al. (2018) and Blanco-Montenegro et al. (2018), suggest that the mantle thermal anomaly beneath the Atlas region seems to be connected with the Canary plume, and that the uplift of the Atlas Mountains (Fig. 1) required mantle upwelling in addition to inversion tectonics and crustal thickening.

As we mentioned in Section 2, petrological and seismic studies propose a thinner lithosphere in the eastern Canary Islands, and for the case of Lanzarote either there is a magmatic layer of underplating up to 8 km beneath the oceanic crust, or Lanzarote is on a region of thinned continental crust. However, the magmatic underplating that changes the strength of the upper crust can create local OTL variations between the islands of about 0.1–0.3 mm. Larger effects of about 1 mm are obtained by changing the elastic properties of layers beneath the Moho boundary up to roughly 220–300 km depths. Park & Rye (2019) note variations in the V_P and V_S velocities of up to 10 per cent at depths of 32 km under hotspots in the Pacific Ocean. The recent DRBD (Debayle et al. 2020) 3-D tomographic model of the upper mantle shows lower shear velocity values for the Canary Islands when compared with those of the earth model STW105 (see Fig. 8).

Bos et al. (2015) demonstrated that changing the shear modulus in the asthenosphere by 10–15 per cent improved the agreement between observed and predicted OTL values in western Europe. We assume that the upwelling under the Canary Islands causes the material to be even less stiff than under western Europe since it is warmer, and thus we lowered the Q values to 75, 50 and 35 in STW105 between the depths 24.4 and 220 km and assumed a Q^α model with α = 0.1. This reduced the V_S velocity by 8, 12 and 16 per cent whereas the shear modulus reduced by 16, 24 and 35 per cent, respectively. As per Bos et al. (2015), the changes in phase-lag are only about 0.2° and will be ignored. Thus, we have constructed new Green's functions based on the STW105 earth model with these specific layers modified (see Fig. 3a) and subsequently computed the respective OTL values. The Green's functions based on Q = 35, 50 and 75 have been included in the Supporting Information. Fig. 9 shows the reductions in the residual differences between observed and predicted OTL for harmonic M₂ using the ocean model FES2014b replaced with CIAM2. For the eastern islands (Lanzarote and Fuerteventura), a reduction of Q values to about 50–35 modifies the misfit between observation and prediction substantially (the average amplitudes vary between 0.6 and 0.8 mm and the phases between 18° and 43° for such reductions of Q). This fact leads us to assume that a decrease of 16 per cent in the S-wave velocities values for the easternmost part of the archipelago results in residuals of about 0.7 mm (Table 4). When we use these modified Green's functions to compute the predicted gravity OTL values discussed in Section 5.1, the misfits with the observed OTL values are still well within the 2σ confidence level of the measurements for Q = 75 and Q = 50 at all stations and on the border for Q = 35 (see Fig. S8). Therefore, upper mantle materials underneath the surface of Lanzarote and Fuerteventura could be less stiff than in the other islands of the archipelago, which is compatible with the low velocity pockets beneath the thin lithosphere of the Canary Islands suggested in Miller et al. (2015), and hence causing the large residuals observed for the semidiurnal band in the eastern islands. This modified mantle, in terms of the elastic properties, underneath the east of the archipelago agrees well with the aforementioned mantle upwelling and is consistent with the fact that the Canary plume is deflected eastward by mantle flow, giving rise to smaller-scale upwelling as suggested by geodynamic models (Duggen et al. 2009; Mériaux et al. 2015; Miller et al. 2015). At the same time, this reduction also causes a slight reduction of the misfit between predicted and observed M₂ vertical OTL displacement for the western islands. This provides an additional argument that a more compliant asthenosphere can produce OTL deformations that are consistent with observations over the whole of the archipelago. In fact, changes in shear modulus, V_S velocity and viscosity could determine the way how the mantle is deformed and, therefore, its dynamics (Irvins & Sammis 1995). A reduction of shear modulus and V_S to adjust the residuals obtained would imply a decrease in the viscosity of the mantle at that depth, which can be explained by an increase in temperature that would lead us to define a thermal anomaly in that zone of the mantle. If we assume that the temperature does not vary and therefore there is no thermal anomaly, the way to lower the viscosity adiabatically is to increase the stress level or deviatoric stress, which implies an increase in the strain rate (e.g. Turcotte & Schubert 2002) and could perhaps represent a change in the deformation regime in the mantle from diffusion to dislocation. However, we can assume (in a simple way) that the asthenosphere in the region has a similar deformation in terms of mechanism and strain rate, and therefore the upwelling of the mantle is due to a thermal anomaly that reduces both viscosity and V_S, which favours the ascent of magmas in that area due to the lower viscosity.

6 CONCLUSIONS

We have used recent global ocean tide models (FES2014b, TPXO9v5a and GOT4.10c), replaced with the CIAM2 regional model, to compare M₂ gravity and GPS-observed vertical displacements with model predictions along the Canary Islands archipelago. We have demonstrated a mostly elastic relation between the weight of the ocean tides and the resulting surface displacement and change in gravity. Also, we note some anelastic effects in the asthenosphere below the archipelago. Indeed, the misfit between predicted and observed OTL for harmonic M₂ is reduced considerably, both for GPS vertical displacement and gravity, when we account for the effects of the anelastic dispersion in the asthenosphere. For gravity, the mean misfit is about 1 nm s⁻² when using the recent seismic models S362ANI and STW105. Furthermore, the gravity residuals support an anelastic behaviour of the upper crust layers in the Canaries, and the use of the CIAM2 regional ocean model to replace the global ones still improves significantly the accuracy of the gravity OTL computation in the Canary Islands region. That is, the discrepancy between the amplitude of the OTL predictions with respect to the observed tidal gravity residuals are reduced from the range −2.2 per cent to 3.5 per cent (without CIAM2) to the range 0.0–0.5 per cent (with CIAM2), at the four observing sites considered.

For M₂ GPS-observed vertical displacement, the average residuals are 0.7–0.8 mm when using FES2014b (or GOT4.10c or TPXO9v5) replaced with CIAM2, and with the STW105 and S362ANI earth models after accounting for anelastic dispersion. However, for some sites on Lanzarote and Fuerteventura, which are located in the easternmost sector of the archipelago and close to continental Africa, the residuals reach up to 2 mm. These high residual values for the eastern islands cannot be explained by errors in the ocean tide models and/or the uncertainty of the GPS observations. However, on modifying the upper mantle layer of 24.4–220 km depth by reducing the value of Q in the asthenosphere, which at tidal periods causes a reduction in the shear modulus of 35 per cent (i.e. assuming less stiff materials in these two layers at tidal periods) through the appropriate anelastic STW105 Green's function, we find a considerable misfit reduction to residuals of about 0.7 mm, similar to that found in the central and westernmost islands of the archipelago. We also demonstrated that this more compliant asthenosphere produces accurate vertical OTL displacement values for the westernmost islands. This is particularly interesting because our results of an elastically more compliant asthenosphere at tidal periods could confirm the existence of the mantle Canary plume, and constrain the most accepted hypothesis of the origin of the Canary Islands that is supported by the unifying model (Anguita & Hernán 2000).

The reduction in the shear modulus by about 35 per cent in the asthenosphere is about three to four times larger than previous findings in other regions of the world, as for instance in western Europe (11 per cent reduction) and East China (8 per cent reduction). This is remarkable because the Canary Islands lie on oceanic lithosphere, whereas the continental lithosphere was involved in the other regions.

ACKNOWLEDGEMENTS

Projects PID2019-104726GB-I00 and CGL2015-63799-P of the Spanish Research Agency has supported this research. Machiel Bos was sponsored by the Portuguese Foundation for Science and Technology in the scope of the project IDL-FCT-UID/GEO/50019/2021 and the ATLAS project (PTDC/CTA-GEF/31272/2017, POCI-01–0145-FEDER-031272 and FEDER-COMPETE/POCI 2020). Nigel Penna was funded through UK Natural and Environment Research Council grant NE/R010234/1. The authors are thankful to the Councils of El Hierro, Gran Canaria and Lanzarote Islands and also to the staff of Izaña Atmospheric Research Center (Tenerife) and the National Park of Timanfaya (Lanzarote) for providing us facilities to carry out the gravimetric observations.Thanks to NASA JPL for the GIPSY software, orbits and clocks." "Authors thank the editor Duncan Agnew, and two anonymous referees for their valuable comments to improve the manuscript.

DATA AVAILABILITY

The GPS data sets are available from stores GRAFCAN (Cartographical Service of Canary Islands Government, https://www.grafcan.es/red-de-estaciones), IGN (Spanish Instituto Geográfico Nacional, https://www.ign.es/web/ign/portal/gds-gnss-datos-rinex), IGS (International GNSS Service, https://igs.org/data-access/#global-dcs). Data from the GPS station CVAN (Gran Canaria) belongs to the Research Group ‘Geodesia’ of the University Complutense of Madrid and is available in the Zenodo repository (https://doi.org/10.5281/zenodo.7324718). The CIAM2 model is available in the Zenodo repository (https://doi.org/10.5281/zenodo.7575945). The gravity data are available from the corresponding author under reasonable request. The tidal analysis software VAV (version 6) and ETERNA (version 3.4) can be downloaded from the International Geodynamics and Earth Tide Service at http://igets.u-strasbg.fr/soft_and_tool.php. The topography and bathymetry data of Fig. 1 are from the GEBCO database (https://www.gebco.net).

CONFLICT OF INTEREST

The authors declare no conflict of interest.

SUPPORTING INFORMATION

Supplementary-Material_Arnoso-etal_GJI-22-0602.R2

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.

References

Author notes