Modelling of crustal anomalies of Lanzarote (Canary Islands) in light of gravity data

Summary

The application of a gravity inversion method enables us to obtain a 3-D density contrast model of the upper crustal anomalies of the volcanic island of Lanzarote (Canary Islands). For this, we use a network of 296 gravity stations distributed over the whole island, and a digital terrain model of about 45000 terrestrial and oceanic data to determine the corresponding terrain correction. A density value of 2480kgm⁻³ is chosen for this correction by means of a new approach. The resulting Bouguer anomaly is analysed by means of a least-squares prediction which gives us a mean level of uncorrelated observational noise of about 1.2 mgal. This anomaly is considered in order to obtain independent information about the inner anomalous mass density distribution by means of a 3-D gravity inversion based on a systematic exploration on a prismatic partition of the subsoil volume, and adopting a priori values of the density contrast (positive and negative) to determine the geometry of the anomalous bodies. The problem of non-uniqueness of the solution is avoided by using a minimization mix condition on the weighted residuals and the weighted whole anomalous mass. The structural solution is finally presented by means of horizontal sections and vertical profiles.

A main intrusive body is located under the central-eastern area and could correspond to a dilated volcanic activity of shield formation. It shows a prismatic form of more than 15km depth, subducted with only the ridges remaining as horst blocks. Moreover, the SW and NE extreme areas of the island show smaller and shallower positive bodies, interpreted as less-developed magmatic intrusions. Conversely, several density lows offer interesting shallow alignments, 45°N (ENE–WSW) and 125°N (WNW–ESE), which could be associated with a fracture system corresponding to structural stress, and also correlate with historic eruptions, such as, for instance, the Timanfaya eruption. The monitoring of several geophysical parameters at two underground geodynamic stations, in the NE zone of the island and Timanfaya, shows characteristic differences between the two zones which confirm crustal anomalies in the second station.

gravity anomalies, inverse problem, Lanzarote island, volcanic structure

Introduction

The Canary Islands (Fig.1) are an old volcanic feature sited on top of Jurassic oceanic crust, located at the edge of the West African Continental Margin. Their origin is still under debate, and authors have proposed different types of genetic models (e.g. Wilson 1973; Schmincke 1982; Holik et al. 1991; Anguita & Hernan 1975; Marinoni & Pasquarè 1994) for the archipelago, where some age determinations (Abdel-Monem et al. 1971) indicate a general E–W age progression. The diverse models apply, for instance, hotspot theory (Wilson 1973), a connection to the Alpine orogeny, which reached its maximum activity in this zone during the Miocene (Anguita & Hernán 1975), and the coexistence of a hotspot and a constraining complex regional structural framework (Schminke 1982). Several geophysical research projects have been conducted in the archipelago (e.g. Bosshard & MacFarlane 1970; Banda et al. 1981; Canales & Dañobeitia 1998; Dañobeitia et al. 1994; Ranero et al. 1997), pointing out clear structural differences among the islands. Among these differences, from recent seismic studies it seems reasonable that the oceanic crust west of El Hierro islands is very shallow, Ranero et al. (1997) deduced a crustal thickness of between 6 and 9km maximum, and beneath the continental shelf there are reported thicknesses of around 20–25km (Hinz et al. 1982; Weigel et al. 1982).

Figure 1

Geographical location of Lanzarote island (29°N, 13.6°W; 635000,3220000 UTM) at the northeastern extreme of the Canarian archipelago (L: Lanzarote, A: Africa). Major volcanic periods in Lanzarote are (Marinoni 1991): (1) the Tableland (series I); (2) the Quaternary volcanism (series II, III); and (3) recent volcanism (series IV). Localities mentioned in the text: F, Famara; M, Mozaga; AJ,Los Ajaches; CV, Cueva de los Verdes; T, Timanfaya; SB,San Bartolomé; PC, Peñas del Chache; AF, Atalaya de Femés; CC, Caldera del Corazonzillo.

The island of Lanzarote is located at the northeastern extreme of the archipelago, in clear alignment with the island of Fuerteventura (Fig.1). The two islands constitute the emergent part of the East Canary Ridge, a NNE–SSW linear volcanic structure located on atypical oceanic crust, at least 11km thick (Banda et al. 1981), between the continental rise and the Canary basin. The East Canary Ridge consists of several uplifted blocks of oceanic basement mantled by a thick pile of volcanic rocks (Grunau et al. 1975), probably overlying an intrusive complex. A sediment cover more than 10km thick is located in the continental rise east of the ridge. Lanzarote is described by a shallow basement, probably about 4–5km thick as deduced from seismic profiles (Banda et al. 1981), formed by a group of sedimentary rocks (quartzite and shales), plutonic rocks (basic and ultrabasic), and subvolcanic rocks (basaltic and trachytic dikes) with an abundance of xenoliths of quartzite and sandstone emitted by its volcanoes (Araña & Carracedo 1978). The lavas are basaltic, and very limited outcrops of massive thrachytes exist in the oldest parts of the island (NW and SE). A general geological study of Lanzarote can be found in Fuster et al. (1968) and Marinoni (1991). Three major volcanic periods are recognized (see Fig.1): (1) the Tableland (series I) basalts dated at 6–12Ma (Miocene–Pliocene); (2) the Quaternary volcanism (series II, III) dated at 1Ma and separated from the former by an erosional interval; and (3) recent volcanism (series IV). Many Quaternary eruptions and several historic ones have been identified. The Lanzarote eruption between 1730 and 1736 was one of the Earth's biggest historical eruptions. The structural evolution of this volcanic island may result from a complex interaction of magmatism with the regional stress field and the local stress field generated during the growth of the island itself.

According to Marinoni & Pascaré (1994), Lanzarote can be divided into three principal morphological blocks (Figs1 and 2). The northernmost part is elongated in a NNE–SSW direction and consists mainly of the Famara massifs. The central part of the island has a characteristic flat truncated conic shape (Mozaga area). The southern block includes the Los Ajaches massifs. These remaining edifices correspond to the dilated eruptive episode of the shield stage (Miocene–Pliocene) with the first subaerial lava flows (Araña & Carracedo 1978). However, large earth movements may have disrupted this central structure until now, and only the ridges remain as horst blocks, the central cone having subsided between them and having been covered by younger volcanic eruptions (MacFarlane & Ridley 1969). Between the pre-erosional and the post-erosional stage, the Tableland succession was deeply eroded and probably tilted about 10° generally to the ESE (Marinoni & Pascaré 1994). After a long period of rest, volcanic activity was renewed (post-erosional stage during the Quaternary) and lasted almost until the present time; the resulting materials and volcanic edifices cover the major part of the island and have been well preserved.

Figure 2

Terrain model of the islands of Lanzarote and La Graciosa. Orthographic SW–NE view. Highest altitudes are about 600m at Peñas del Chache (NE Lanzarote). In this figure two ancient and eroded edifices, Los Ajaches in the south and Famara in the north, are observed. One third central massif appears more eroded and subducted with respect to the former. In this central zone many well-preserved volcanic cones can be observed.

So, the present structural architecture is the result of a complex magmatic and tectonic evolution characterized by variation in the stress field from the Miocene to the present. According to Armenti et al. (1989), the present structure is largely controlled by post-Miocene displacements which divide the island into five NW–SE elongated blocks, represented by the three structural highs (Famara, Mozaga an Los Ajaches) and two grabens. These blocks were subsequently dissected by NNE–SSW and NE–SW faults. The greatest principal stress shows an ENE–WSW direction, which argues for a transpressive regime. The NE–SW-striking faults were generated by extensional tectonics with a minor transcurrence component. This resulted in a larger down-throwing of the northwestern sector relative to the southeastern one. A third narrow and elongated central graben, marked by the Timanfaya volcanic alignment, could be delineated between these two sectors.

Most of the geophysical works about the Canary Islands correspond to extended areas and marine data. Here we look for a local study corresponding to the upper crust of Lanzarote. To obtain a better understanding of the inner structure and evolution of Lanzarote, several geophysical and geodynamical studies were carried out on the island by the Instituto de Astronomı´a y Geodesia (IAG, Madrid). Two underground geodynamic stations monitoring several geophysical parameters were set up (Vieira et al. 1991; Arnoso et al. 2000). The station Cueva de los Verdes is located inside a lava tunnel of the La Corona volcano (Fig.1) at the NE end of the island. The second station is located above the geothermal anomaly area in the Timanfaya National Park (Fig.1). Shukowsky & Mantovani (1999) suggest that lateral variations in the scales of tectonic units could be inferred from tidal records. Continuous gravity tide observations have been carried out in these two laboratory locations for several years, and have been used to obtain the respective local tidal gravity models (Table1), which help in studying the anomalous crustal structure of the island in these two zones (Arnoso et al. 2000).

Table 1

Observed amplitudes A and local phases α for the main tidal waves O₁ and M₂ at stations Cueva de los Verdes (coordinates: 29°.16 N, 13°.44 W, 37 m) and Timanfaya National Park (coordinates: 28°.99 N, 13°.75 W, 381 m). Amplitudes are given in µgal (10⁻⁸ ms⁻²) and phases in degrees.

The IAG has also carried out several gravity studies on the island (e.g. Sevilla & Parra 1975). In this sense, gravity modelling plays an important role in studies of volcanic structures. Gravimetric methods have been employed by several authors (Malengreau et al. 1999; Rousset et al. 1989; Loddo et al. 1989; Camacho et al. 1997; Rymer & Brown 1986), and prove to be useful in the study of the origin, structure and activity of volcanic areas. The work of MacFarlane & Ridley (1969) corresponds to a smaller land gravity data set for Lanzarote. Gravity investigations do have, however, certain well-known limitations, the two most prominent perhaps being the deep resolution and the non-uniqueness of the solution. These disadvantages can be reduced be means of appropriate tools (geological information or statistical techniques).

In 1988 a detailed gravity survey was carried out on Lanzarote island, the analysis of the data and preliminary structural model being presented in successive papers (Camacho & Vieira 1991; Camacho et al. 1992). The aim of the present paper is to obtain an improved interpretation of this gravity data by means of a 3-D inversion approach, to determine the geometry of the anomalous bodies corresponding to prescribed densities. This inverse method involves the systematic exploration of the model possibilities, admitting simultaneously positive and negative contrast (Camacho et al. 2000).

Gravimetric and topographic data

The gravity data consist of 296 stations covering the islands of Lanzarote and La Graciosa with a nearly homogeneous distribution and a data separation of 1.5km (Fig.3). A LaCoste-Romberg gravimeter with digital electronic reading was used in the observation work, and a total of 437 observations were registered.

Figure 3

Distribution of gravity stations and map of the Bouguer gravimetric anomaly. The value used for the density of terrain correction is 2480kgm⁻³. Contour interval are 5mgal. UTM coordinates in metres.

First, the gravity data were corrected taking into account the tidal model calculated using the records of the tidal gravity stations on the island (Arnoso 1996). Then, using the redundant observations, a global least-squares adjustment was made, taking as unknowns linear drifts, possible record jumps (as a result of mechanic or thermal shocks), and the (incremental) gravity values for each station. This adjustment produced residues with a standard deviation of 0.136mgal (1mgal=10⁻³ms⁻²). Finally, the adjusted gravity data were referred to IGNS71. The station coordinates were determined by locating the points on detailed 1:5000 charts.

To calculate the further terrain correction of the gravity data, a digital terrain model was obtained for Lanzarote, La Graciosa and surrounding oceanic areas by means of a dense regular digitalization of 1:25000 topographic charts. Fig.2 shows an orthographic view of the resulting terrain model (45000 data points). The highest altitudes are about 600m. The figure shows two characteristic terrain elements. First, there are the two extreme areas of strongest relief: the Massif of Famara (northern) and the Massif of Los Ajaches (southern). These remaining edifices correspond to a dilated eruptive episode and posterior processes of large earth movements, subsidence, volcanic eruptions and an erosional stage. Second, there are the volcanic cones in the central area of Lanzarote and in La Graciosa. These preserved edifices, corresponding to the post-erosional stage, follow characteristic alignments.

Because of the proximity of Fuerteventura island, its contribution to the terrain correction was recalculated taking into account the terrain model of this island presented by Montesinos (1999).

Gravity anomaly

Using the 1980 normal gravity formula, and determining the terrain correction (extended to 45km, according to the regularity of the effect beyond this distance), the refined Bouguer anomaly Δg was determined (Fig.3). In choosing the density value for the terrain corrections we propose a new approach.

Several methods to calculate the terrain density, such as that of Nettleton (1939) which minimizes the correlation between the Bouguer anomaly and the topography, are efficient if the density of the formations is not related to the subsoil anomalous structure. In volcanic terrain, however, this correlation is generally predominant. The highest areas, for example, often correspond to the eruptive centres, which involve density anomalies.

An inadequate terrain density would assign an excess of topographical anomalous masses as sources of the anomaly field, producing some fictitious local gravity anomalies. Starting from the gravity data themselves, we look for a suitable terrain density that minimizes the gravity anomalies of topographic origin, keeping the anomalies due to the deep bodies. A tool to detect the presence of gravity anomalies of topographical origin can be obtained by means of the anomalous masses model resulting from the inversion process. In fact, an inadequate terrain density will produce a fictitious increase of the local anomalies. It will cause an excess of (shallow) anomalous masses as sources of the anomaly field. Therefore, the proposed criterion is to choose the density that produces a minimum value of the whole anomalous mass of the model that results from applying the inversion approach.

To apply this criterion, however, we must take into account that, in the inversion process (see next section), the whole resulting magnitude of the anomalous mass also depends on the fit level (mean resulting residues) selected for the gravity anomalies (which, in turn, is regulated by means of a factor in the inversion process). Thus, to reach this desirable minimum of mass, we test the gravity inversion results corresponding to several terrain densities and to several possible fit levels (mean residual level).

This method was applied to the gravity anomalies of Lanzarote, resulting in a mean value of 2480kgm⁻³ for the terrain mass density. Fig.4 shows the results of various inversion tests (according to the inversion method detailed below) corresponding to several possible terrain correction densities and to different allowed mean residues. The chosen value is an intermediate density somewhat greater than 2300kgm⁻³, suggested by MacFarlane & Ridley (1969) as a realistic value for the upper layer of the crustal structure from seismic refraction data, but smaller than the 2560kgm⁻³ that can be obtained by Nettleton's method (null correlation altitude-anomaly).

Figure 4

Relation between resulting whole anomalous mass (×10⁹ kg) and various allowed mean residues (standard deviation, mgal) for several inversion tests corresponding to several possible terrain correction densities. The value of 2480kgm⁻³ provides the minimum whole anomalous mass and is therefore selected for the terrain correction.

With the adopted terrain density, the terrain corrections range between 2 and 12mgal with a mean value of 3.3mgal, and the resulting Bouguer anomaly ranges between 152 and 223mgal (Fig.3).

A study of the Bouguer anomaly by means of covariance analysis of the data and a further least-squares prediction (e.g. Montesinos et al. 1999) enable us to filter the data and obtain a value 1.2mgal as the standard deviation of the non-correlated noise (Fig.5). The value is conditioned not only by the accuracy of the (gravimetric and topographic) data, but also by the distance (about 1.5km) between stations and the horizontal gravity gradients (about 3mgal km⁻¹) present in the area. This method enables us to obtain a suitable covariance matrix, corresponding to the inaccuracies for the gravity data, that is useful in the following inversion process.

Figure 5

Histogram of the residues of the filtering process. This non-correlated noise of the data follows a normal distribution of mean 0mgal and standard deviation 1.19mgal.

Gravity inversion method

While limited by the validity of the adopted hypothesis, the gravity inversion methods that seek to determine the geometrical properties of anomalous bodies with fixed density contrast correspond to a non-linear environment and offer interesting results. Usually these methods (linearized approaches, iterative gradient methods, etc.) start from good initial solutions and improve iteratively (e.g. Barbosa et al. 1997; Enmark 1981). General complications come from adopting a 3-D environment and when positive and negative anomalous contrasts are simultaneously accepted. In addition, the absence of a good initial model combined with data inaccuracies and uncertainties concerning regional trend introduce more ambiguity into the inverse problem. For the fully non-linear treatment, the methods of exploration of the space model are often the best option (Tarantola 1987). This exploration process can be conducted randomly (Silva & Hohmann 1983) or by means of a systematic approach, according to the size of the model space to explore. Here we consider an inversion method, a modification of the one described in Camacho et al. (2000), which can be included in this last group of systematic exploration.

For this non-linear method, the usual strong hypothesis is adopted: the subsurface anomalous structure is characterized by prescribed mass density contrasts. Moreover, positive and negative density contrasts are simultaneously accepted. Therefore, the problem consists of determining the geometry of the anomalous volumes corresponding to those prescribed density contrasts, and responsible for the observed gravity anomaly. To describe that geometry, the subsurface volume is divided into a fixed discrete 3-D partition of parallelepiped cells, and then the anomalous volumes are constructed and described as aggregations of elements filled with the prescribed density contrasts.

The main problem for the applicability of the exploratory methods is the large size of the space model to be explored. In our case, a direct application (test across every possible geometrical configurations) of the systematic exploratory search would be very tedious. This is solved by using a ‘growth’ process or expansion approach to construct the anomalous bodies. In this case, the exploration of every model possibility is substituted by exploration of several possibilities of growth (cell by cell) for each step of the expansion of the anomalous bodies. This is a more reasonable goal. Thus, step by step, the prismatic cells, each with prescribed densities, are systematically tested, and then the best options are adopted and added to the growth approach for the anomalous bodies construction. Simultaneously, a regional trend g_reg(x,y) of the gravity data can be adjusted.

For the kth step of the growth process, the residues v_i, i=1,…,N, for the N gravity stations can be calculated as

 

1

where Δg_i is the observed anomaly; A_ij is the vertical attraction of the jth cell on the ith gravity station for unit density; Δρ_j is the density contrast value assigned to the jth cell (among the possible prescribed values); J_k is the set of indexes corresponding to the cells filled until this kth step and then constituting the growing anomalous body; x_i,y_i are the station coordinates; and f≥1 is an adjusted scale factor for the fit.

If only a minimization criterion for the residuals v is used, the acceptance of positive and negative values for the prescribed density contrasts and the inclusion of the trend unknowns give a non-uniqueness problem. To remedy this, an additional condition of minimization of the model variation can be adopted. Thus, the solution is obtained by a mixed condition formed by the gravity fit with the l₂ norm and the whole anomalous mass quantity, using the parameter λ for a suitable balance:

 

2

where v=(v₁,…,v_N)^T are the gravity residues given by (1) for the N stations, m=(Δρ₁, …,Δρ_NB) are the anomalous densities for the NB cells of the subsoil partition, λ is a positive factor, empirically fixed, for balance between model fitness and model smoothness, and f≥1 is the adjustable scale factor which keeps the model gravity nearly proportional to the observed value. Q_D is the covariance matrix (usually a diagonal matrix) corresponding to the estimated (Gaussian) inaccuracies of the gravity data, and Q_M is a covariance matrix corresponding to the supposed determinability of the model parameters m. We take, as covariance matrix Q_M, a diagonal normalizing matrix of non-null elements the same as the diagonal elements of A^TQ_D⁻¹A. The λ parameter governs the application of the minimization conditions regarding the balance between total anomalous mass and residual values (Fig.4). For low λ-values a better fit is obtained, but the anomalous mass may increase excessively and adopt too deep a position. For high λ-values the adjusted model can be too poor and shallow. The λ-value is chosen so that the standard deviation of the inversion residues is similar to the data noise level (1.2mgal for the data of Lanzarote) determined by the previous covariance analysis.

For each step of the anomalous growth, a systematic exploration of the expansion possibilities is tested with respect to the value given by (2), determining also the corresponding f-value and the parameters of a regional component adjusted as a linear trend. This expansion process stops when the scale factor f is close to one. Finally, the solution appears as a 3-D distribution of prismatic cells that have been assigned some of the prescribed contrast densities. Moreover, a regional trend is also obtained, complementing this anomalous mass distribution.

This inversion method, detailed in Camacho et al. (2000), requires, as usual for many non-linear methods, a suitable choice of the a priori density contrasts (positive and negative in our case). A too-high density contrast would fit the main anomaly component but would produce a rather ‘skeleton’ model, perhaps unable to represent very slight anomalous structures. A too-low density contrast would fit the light structures, but would produce a rather ‘inflated’ model, perhaps unable (taking into account the volume limitations) to fit the main components, and giving rise to fictitious peripheral bodies. To remedy this, and if geological knowledge of the true local density contrasts is not available, we propose an improvement of our inversion method by means of a theoretical approach.

Instead of invariable density contrasts along the whole model formation, we propose now to adopt changing contrasts. Starting from fixed maximum values (let us name them R_P positive and R_N negative), the prescribed density contrasts evolve during the growth process according to a simple law. By means of empirical tests, we have selected for each growth step the density contrasts Δρ_P, and Δρ_N, positive and negative, given by

 

3

where f≥ 1 is the scale factor corresponding to the step of the fit process, and τ is a fixed factor corresponding to the desired variability of the density values. A high τ-value will produce an anomalous model with contrast values mostly close to the maximum values (and thus with a sharp geometry). A low τ-value will produce a model with more variable or decreasing contrasts (and thus with a more diluted geometry). In this form, the high density contrasts (close to the maximum values R_P and R_N) will be employed to fit the main anomaly components, while smaller density contrasts fit the smaller and local anomalies. The resulting anomalous model is richer (and perhaps somewhat more diffuse).

We must keep in mind that the gravity inversion is nearly insensible to horizontal mass stratification, and therefore the obtained model must be considered independent of any additional stratification. In spite of these limitations, this method provides interesting gravimetric information. The method has the following advantages:

• 1

  it adopts a 3-D environment;

• 2

  non-gridded, non-planar and inaccurate data are accepted;

• 3

  a simple regional trend can be simultaneously determined;

• 4

  adjustment of an indefinite number of anomalous bodies is possible;

• 5

  non-regular subsoil partition can be considered (e.g. with deeper blocks bigger than shallow blocks);

• 6

  if previous evaluated models exist, they can be incorporated;

• 7

  a variable density contrast can be adopted in the inversion process; and, above all

• 8

  positive and negative density contrasts are simultaneously accepted (see Camacho et al. 2001 for software development).

3 Adjusted -d model of anomalous density

Taking into account that the distance between contiguous stations is about 1500m and that the mean diameter of the survey is about 30000m (and the computer limitations), we adopted a 3-D partition of the local subsoil into 16800 rectangular prisms, with sides ranging from 700m for the shallow elements (1000m depth) to 1300m for the deepest blocks (at 18000m depth). According to the suggested method, we selected, for the whole subsoil volume, the possible density contrasts by means of the extreme values, R_N=−440kgm^–3 and R_P=+440kgm^–3, and an exponential variation law (with τ=4) of the contrasts. These values have been selected empirically to obtain connected and well-developed volumes and to produce an inversion model with a negative mean density contrast of −346kgcm⁻³ and a positive one of 341kgcm⁻³, which could correspond to the usual contrasts among structural elements of the Canarian volcanism (these values are close to the contrasts that occur when considering, for instance, density values of 2250kgm^–3 for sedimentary material, 2700kgm^–3 for crustal material, and 3100kgm^–3 for mantle-type material) (Bosshard & MacFarlane 1970). When different density contrasts are used, the resulting solutions are rather similar, although for higher contrasts, smaller volumes would be obtained.

Given this partition of the subsoil, and considering the above density contrasts, we applied our inversion method and obtained a model with the geometry of the anomalous bodies (for the adopted densities) and a simple regional component of the gravity anomaly field by means of a polynomial surface of degree one. This component reflects the general tendency of the anomaly and is strongly conditioned by the two main anomaly highs located in the east and northeast of the island.

The adjusted local model of density contrasts is shown by several horizontal sections and vertical profiles in Figs6 and 7, from which some interesting features can be highlighted.

Figure 6

Adjusted model of anomalous density contrast obtained by means of 3-D gravity inversion. Several horizontal sections, from −2000to −14000m depth, are shown.

Figure 7

Adjusted model of anomalous density contrast obtained by means of 3-D gravity inversion. Several vertical profiles across the main structures (according to the lines shown in Fig.6d) are shown.

The main structural element in the adjusted model is a positive contrast mass (‘A’ in Fig.8) which emerges close to the centre of the island, in the area of San Bartolomé (see Fig.1). Its lateral size is about 13km (Fig.6), and, according to the adjusted model, its anomalous mass is about 582×10¹² kg. The top of this body, located in the shallower sections of the model, presents a double summit (‘A1’ and ‘A2’ in Fig.9) with both summits elongated along 45°N and with a mutual strike of 125°N. Its bottom appears in the adjusted model with a depth greater than 15km (Figs6 and 7) (not far from the supposed depth of the Moho). The maximum depth of this body is, however, somewhat doubtful, because its determination interferes with the regional trend adjustment. While for the deep sections (under 7km) this main body appears rather rounded in the model, for mean depths (1–6km), where the gravity survey is more sensitive, the model shows a clear geometry with a nearly square form according to directions 45°N (NNE–SSW) and 125°N (NNW–SSE). This important anomalous mass with its nearly prismatic shape can be interpreted as an intrusive body. According to MacFarlane & Ridley (1969), this body corresponds to a major igneous centre within the crust; it was important at least during the early subaerial growth of the island, and may have formed initially at the intersection of fundamental NNE- and SW-trending fault systems. However, while large earth movements and an erosional process may have disrupted this structure until now, only the ridges remain as horst blocks, the central cone having subsided between them and having been covered by younger volcanic eruptions. This hypothesis is in agreement with a possible interpretation of our model taking into account the irregular form (with negative and positive contrasts) in the shallower sections, associated with recent volcanism and the presence of the two structures A1 and A2. Marinoni & Pasquare (1994) relate this intrusive body to surface outcrops of the Shield phase. An alternative explanation of the very superficial low between A1 and A2 could be a simple erosion process.

Figure 8

Horizontal section of the anomalous density contrast model corresponding to a model mean depth of 4000m. The fairly regular geometry of the positive density contrast bodies is highlighted. These structures are interpreted as a system of intrusive bodies with SWS–ENE trend. Subsidiary bodies thus appear as drained from the main sources (the central body A and the smaller structure B).

Figure 9

Shallow horizontal section of the anomalous density contrast model corresponding to a depth of 2000m. A pattern of alignments of negative density contrast zones, with main directions ENE–WSW and WNW–ESE, is identified and associated with fracture areas and recent volcanism.

A subsidiary intrusive body is located from 1 to 6km depth close to the main body (body B in Fig.8), with an anomalous mass of about 33×10¹²kg. This structure suggests a much smaller accumulation of magma in the crust, drained from the central reservoir (A in Fig.8) along a SSW–NNE-trending fissure system indicated as the alignment of the two bodies. The system follows the WSW–ENE direction of the principal stress in the island (Armenti et al. 1989). The body is situated to the south of the Peñas del Chache, in the Macizo de Famara (see Figs1 and 2), the highest point of the island and one of the oldest areas.

Two positive bodies are also noticeable: a peripheral intrusive body (C in Fig.8), located SW of Lanzarote and elongated along 126°N; and a third intrusive body (D in Fig.8) located at the NE extreme of Lanzarote and elongated along 133°N. This last body is too peripheral for its limits to be determined, although the adjusted model suggests a magnitude smaller than that of the main body. The gravity data observed in La Graciosa Island do allow us, however, to define the geometry of this body on its SW side.

Body C lies west of Atalaya de Femés volcano (see Fig.1), which was constructed in the post-erosional stage (Marinoni & Pasquaré 1994). It extends between the shallower sections of the model and 7km depth with an anomalous mass of 48× 10¹²kg. A more superficial and smaller body appears aligned along 43°N (E in Fig.8), close to several small outcrops of Series I (Coello et al. 1992).

The existence of these two less developed and shallower positive bodies, C and D, suggests that the accumulation of mantle-type material took place close to the surface. Between these not very massive SW and NE areas, the central body, A, suffered a further subsidence process during growth (Watts 1994), producing an extensive body with deep roots as is observed in the deeper sections. According to this interpretation, MacFarlane & Ridley (1969) note that the topography and geology of the island suggest that the central part has been down-faulted between the uplands in the north and south, forming a graben structure.

All these positive bodies present a nearly regular prismatic shape (within the intrinsic smoothness of shape for the deep zones in the inversion of potential fields), with borders along 45°N and 125°N. Moreover, the adjusted model shows others characteristic alignments. The main body, A, shows a like-alignment with the bodies B and D, with strike about 23°N. Meanwhile, the bodies in the south of the island (C and E) lie, with similar strike, along a line 15km displaced to the west, and they could be connected by alignment with the crustal structures identified in the north and central area of the neighbouring island Fuerteventura by means of a similar gravimetric interpretation (Montesinos 1999). So, MacFarlane & Ridley (1969) also suggest lateral movements as inferred by the NNE-trending scrap visible in the northern uplands that appears to have been displaced 15–20km westwards in the south of the island.

Two main kinds of low-density areas are shown in the adjusted model: extensive, peripheral and deep minimum structures that limit the E and SE borders of the island; and the internal minimum areas which form elongated trenches.

The first group, extrapolated in the model, and mostly on the west side of the island, may be associated with oceanic sediments coming from the nearby African margin that infill the moat as a result of the subsidence of the main intrusive body A. We observe (Figs6 and 7) that these structures appear from the surface down to a depth of about 7km. However, taking into account the extrapolated nature of these lows, reliable geometric properties cannot be deduced.

The second group is composed of thin elongated structures that run across the centre and western area of the island with depths ranging from 1 to 4km. These trenches follow two main directions, 45°N(ENE–WSW) and 125°N(WNW–ESE), across the island (Fig.9). The Timanfaya fissural eruption occurred (1730–1760) within this minimum-density area; its fractures and volcanic alignment followed two main directions, E–W and ENE–WSW, each one defined by smaller alignments of superposed vents (Marinoni & Pasquaré 1994). Connected with this trench pattern, we observe a main minimum of greater depth (Fig.8) under the cone of Caldera del Corazoncillo (see Fig.1), in the zone of Timanfaya, at the intersection of two alignments of minima according to the main directions described (Fig.9).

Other peripheral trenches of low density with local directions are also identified in the model; for instance, with a direction of nearly N–S in Famara, and a direction of nearly E–W in the northwestern part of the island (Figs6 and 9). In La Graciosa, another small local minimum is identified in the eastern part of the island. We consider all these elongated low-density structures to be connected with fractured areas and recent volcanism.

The obtained model is in agreement with the hypothesis that the central area of Lanzarote is a subducted structure, crossed by a system of shallow fractures, indicating a crustal difference with respect to peripheral formations at Famara and Los Ajaches (e.g. Armenti et al. 1989). The NW zone of this central area, and especially Timanfaya, appears to be subjected to the main fractures deduced from the adjusted model. These crustal differences have also been observed by means of several techniques that highlight discrepancies between results obtained in Timanfaya and other results in the zone of Famara (Laboratory of Cueva de los Verdes). For instance, the gravity tide anomalies computed after correcting the ocean loading and attraction effects (Arnoso et al. 2000) show clear differences in magnitude at the two stations, being larger for M₂ and O₁ at Timanfaya than at Cueva de los Verdes station. Moreover, the phases of the residuals for O₁ and M₂ are more consistent with a body tide effect than with an ocean loading effect, and the sign of the O₁ and M₂ residuals at Timanfaya station are negative, which indicates a response of a porous or fractured local upper crust subject to the influence of tidal strain (Arnoso et al. 2001). Dı´ez Gil et al. (1986), by means of thermal anomalies, suggest the existence of a magmatic intrusion, of radius 200±100m, located at a depth of 4±1km, with a convective system in which the energy is transported by fractures which continue to depth according to electromagnetic studies. These results are compatible with air model and with the existence of bodies of negative density contrast such as those indicated.

Finally, Fig.10 shows the distribution of residues of the inversion adjustment, which have a standard deviation of 1.15mgal, in accordance with the data accuracy.

Figure 10

Distribution of residues corresponding to the inversion process. A standard deviation of 1.15mgal, close to the value obtained for the data noise, is presented.

Conclusions

The application of a new version of the gravity inversion method proposed in Camacho et al. (2000) on the gravity data observed in Lanzarote and La Graciosa enables us to obtain a model of crustal anomalous mass for these islands. The gravity values considered have been analysed, resulting in a level of data accuracy of 1.2mgal (non-correlated noise). To obtain the Bouguer anomaly we use a density for terrain correction of 2480kgm⁻³. This value is calculated by means of a new approach based on the whole anomalous mass resulting for the inversion model. The obtained Bouguer anomaly ranges between 152 and 223mgal, and is used to apply the modified inversion method. The actual improvement of the method consists of permitting the variability of the density contrasts that appear in the model. The adjusted model for Lanzarote is constructed by limiting the possible density contrasts to between −440 and 440kgm⁻³, resulting in a geometrical configuration of high- and low-density areas. The results show a main intrusive body as a result of the accumulation of mantle-type material and further subsidence. Other smaller and shallow intrusive bodies are detected. These positive bodies can be related to the shield stage of the island formation. Their geometry is characterized by regular shapes with clear strikes (45°N and 125°N). This regularity suggests a tectonic component in the origin and evolution of these structures. The adjusted low-density bodies comprise some peripheral structures and several linear trenches which we associate with a regular fracture system. Some of these crustal anomalies have been confirmed using results obtained in two geodynamic laboratories (in Timanfaya and Cueva de los Verdes). The gravity tide anomalies show a larger magnitude at Timanfaya than at Cueva de los Verdes station, and the results at Timanfaya station indicate the existence of fractured and porous upper crust.

Acknowledgments

This work is part of the projects PB88-0022, PB91-0150 and AMB97-0706 conducted by the Instituto de Astronomia y Geodesia and developed in the Canary Islands. Additional support was obtained from the Cabildo Insular de Lanzarote. J. J. Dañobeitia provided helpful comments on the manuscript.