Abstract

We propose a method to evaluate the stress generated at the local scale by the spatial variations of the gravitational potential energy (GPE), which is related to inhomogeneous topography and mass distribution in the lithosphere. We show that it is possible to infer these local stress sources from the second spatial derivatives of a geoid height grid, used as a proxy of the GPE.

The coherence of the method is validated on a passive margin, the Bay of Biscay. The result is that expected in such a geological configuration, with extensional local stress sources with the maximum horizontal principal stress parallel to the margin and compressive sources with the maximum horizontal principal stress perpendicular to the margin in the continental and oceanic lithosphere, respectively.

We apply the method to Western Europe in order to provide a better understanding of the complex spatial variation of the present-day tectonic activity. Our results indicate a stress pattern from the local sources dominated by short-space-wavelength (of the order of a few tens of kilometers) variations in the tectonic style and in the direction of the maximal horizontal principal stress σH. A comparison of the σH orientations and tectonic style from the local sources with the ones of the World Stress Map (WSM) data set indicates that the local stress sources can be representative of the deviatoric stress state in some regions. Our results explain 71% of the faulting styles for the earthquake fault-plane solutions in the WSM, which is better than the classical compressive NW-SE stress field model. In the central part of the Pyrenees, the agreement between earthquake fault-slip directions and the direction of shear stress from the local sources acting on the associated fault planes is compatible with the extensional stress field evidenced by recent investigations.

INTRODUCTION

The main source of data worldwide about the state of stress within the lithosphere is the World Stress Map database (WSM), a global compilation of local stress indicators (Zoback et al., 1989; Zoback, 1992; Heidbach et al., 2010). From the first analysis of the WSM database release, Zoback (1992) suggested the existence of long-spatial-wavelength, first-order compression stress patterns in plate interiors that are related to plate driving forces, mainly ridge push and continental collision. In addition, buoyancy forces generated by large-scale lateral variations of the density in the lithosphere and surface topography are capable of modifying this first-order stress pattern (Artyushkov, 1973; Fleitout and Froidevaux, 1982, 1983). Zoback (1992) also mentioned that regional or second-order stress fields resulting from flexural stresses and local crustal strength can result in the rotation of regional stress orientations. Heidbach et al. (2007, 2010) suggested that, even if the global stress pattern is primarily controlled by plate-boundary forces, local and regional stress sources can be of the same order of magnitude as the global one, particularly in Europe, where short-wavelength (smaller than 200 km) stress patterns are frequently observed. Despite the richness of the WSM, information remains limited to specific regions, and fundamental questions remain unsolved, not the least of which is to explain the seismic activity in plate interiors.

In parallel to the information provided by the WSM database, the state of stress in the lithosphere is also investigated by modeling. In some studies (Flesch et al., 2007; Ghosh et al., 2009; Naliboff et al., 2012), the stress field in the lithosphere is evaluated by using the thin sheet approximation for the lithosphere and the associated horizontal force balance equations (England and McKenzie, 1982). The horizontal spatial derivatives of the gravitational potential energy (GPE), equivalent to the depth-integrated lithostatic pressure in the lithosphere, which are the first type of body-forces–like terms of these equations, are evaluated either by using lithospheric thicknesses and densities obtained by independent geophysical determinations, or by using the geoid as a proxy for the GPE. Most of these recent studies are large-scale investigations (Flesch et al., 2007; Ghosh et al., 2009; Naliboff et al., 2012), and few authors (Pascal and Cloetingh, 2009) have conducted more regional studies. The second type of body-forces term in the force balance equations is the horizontal tractions at the base of the lithosphere arising from deep mantle convection, at the origin of the dynamic topography. Their contribution to the stress state can be as important as the one of the GPE in areas of continental deformation (Ghosh et al., 2008; Faccenna and Becker, 2010).

In Western Europe, the stress state is classically attributed to a compressive NW-SE far-field stress due to plate-boundary processes: Europe-Africa convergence and ridge push from the Mid-Atlantic Ridge (Grünthal and Stromeyer, 1992). However, from the 958 data points with quality A to D included in the WSM, only 532 present σH directions ranging between 105° and 165° (135° ± 30°), meaning that the majority, 55.5%, of σH is oriented around the NW-SE direction, with a variance of 30° that is equivalent to a C-quality uncertainty; more details on the ranking of the WSM data quality can be found in Zoback and Zoback (1989) and Sperner et al. (2003). The WSM data set includes 542 focal mechanisms; 40% of them present an extensive style that cannot be explained by the compressive far-field stress. Hence, as suggested by Heidbach et al. (2007, 2010), local stresses can control the stress pattern at different spatial scales in Western Europe. This is also evidenced by the small-scale spatial variability of earthquake fault-plane solutions (Delouis et al., 1993; Camelbeeck and van Eck, 1994; Kastrup et al., 2004; Mazabraud et al., 2005).

This is why we investigated the ways in which the short-wavelength spatial variations of the GPE in Western Europe contribute locally to the stress state in the lithosphere. For that purpose, we used the geoid in place of the GPE evaluated from topography and crustal model because, in addition to the fact that the geoid is obtained from observation rather than modeling, we assume that the spatial resolution of the geoid, around 5 km in Western Europe, is better than the modeled GPE and that, at the short length scales of our study, the long-wavelength contributions resulting from sublithospheric buoyancy forces are unresolvable and can therefore be ignored.

By only considering the stress contribution from the GPE, we ignored the part of the stress state resulting from the tractions at the base of the lithosphere caused by the mantle flow. On the one hand, Western Europe is a region where global positioning system (GPS) measurements, as for a large part of Europe north of the Alps and Carpathians, show that it behaves rigidly and defines a stable reference frame (Nocquet and Calais, 2004; Nocquet et al., 2005). This indicates that mantle flow has little influence on the local stress state, despite the fact that predicted dynamic topography maps (Faccenna and Becker, 2010) suggest significant mantle flow in the Pyrenees and the Bay of Biscay. On the other hand, as already mentioned, part of the WSM local stress observations in Western Europe cannot be explained only by large-wavelength stresses. Therefore, an evaluation of the local sources of stress will allow a comparison with the WSM observations and discussion of their relative importance in the light of the different stress contributions in the lithosphere, including the one from the mantle flow.

The developed method, presented in the following two sections, evaluates the integrated stress generated in each column of lithosphere on a grid of 5 km × 5 km from the second spatial derivative of the geoid height. The method is tested in the subsequent section on the case of a passive margin, the Bay of Biscay, a geological context for which the stress distribution caused by the contact between oceanic and continental lithosphere is well understood. We evaluate these local sources of stress in Western Europe and in the Pyrenees (in the last two sections) and their relative importance to the local stress state by comparison with the World Stress Map database for Western Europe and earthquake fault-plane solutions for the Pyrenees. Figure 1 shows the three investigated areas.

STRESSES IN THE LITHOSPHERE AND GEOID HEIGHT GRADIENT

If the shear stresses at the top and at the bottom of the lithosphere are negligible, the vertical stress σ(z)zz at depth z in the lithosphere is a principal stress, and σxz and σyz are negligible. According to Ghosh et al. (2009), tractions at the base of the lithosphere arising from mantle convection are small enough that vertical may still be considered as a principal direction. The associated principal stress value is equal to the weight per unit area of the overlying rock (Artyushkov, 1973; England and McKenzie, 1982, 1983; England and Jackson, 1989): 
graphic
where ρ(z) is the density at depth z, and h is the surface elevation.
The averaged value of the lithostatic pressure is given by: 
graphic
where L is the depth of the base of the lithosphere, and GPE is the gravitational potential energy per unit area (Ghosh et al., 2009).

Both lateral differences of density, as well as of surface topography and its compensation at depth in the lithosphere, produce lateral differences in σzz and of the buoyancy force. Those lateral differences give rise to horizontal stresses in the lithosphere acting from the region where the integral of σzz is high toward the region where the integral of σzz is low.

The deviatoric stresses τ in the lithosphere read as: 
graphic
where p = 1/3(σxx+ σyy+ σzz) is the pressure, σxx and σyy are the normal stresses on surface perpendicular to the horizontal reference axes, and σxy is the horizontal shear stress.
We define the average deviatoric stresses in the lithosphere as forumla and we have the relationship: 
graphic
In the thin sheet approximation (England and McKenzie, 1982, 1983), the horizontal force balance requires the equality of the buoyancy forces per unit area with the vertically averaged deviatoric stresses forumla which gives: 
graphic
 
graphic
In the one-dimensional (1-D) case of a vertical plane in the x and z directions, the force balance equations are reduced to one equation with no shear stress component: 
graphic
In this case, from Equation 3, we obtain forumla and the integration of Equation 6 gives: 
graphic

The horizontal stress in the lithosphere column at x is proportional to the difference of the column gravitational potential energy per unit area (GPE), where GPE0 is the value of a reference column in which the stress state is assumed to be zero forumla in the absence of far-field forces (Houseman and England, 1986). The deviatoric stress acting in the column is a tensile pressure forumla and forumla, or a compression pressure forumla and forumla, if the difference GPE(x) – GPE0 is positive or negative, respectively.

Turcotte and Schubert (2002) provided an expression for the horizontal forces per unit length F in an isostatically compensated lithosphere due to the lateral differences of density, as well as of surface topography and its compensation at depth as a function of the surface geoid anomaly ΔN: 
graphic
Hence, the difference of pressure from a lithospheric column at (x, y) to the neighboring columns of lithosphere at (x′, y′) is given by Equation 8, where forumla with forumla is the unitary vector in the direction of (x′, y′). The horizontal body force density forumla inside the column is proportional to the components of the geoid gradient 
graphic

The two components of the depth-integrated horizontal force in the column are given by FxΔx and FyΔy. They are expressed in N/m.

By interchanging the differential and integral operations in Equation 8, we obtain a relationship between the horizontal forces as deduced from the geoid anomaly between two surface points and the difference of the GPE at those points (Artyushkov, 1973; Jones et al., 1996). 
graphic

To be valid, the limits of the integral, –h and L, should remain unchanged. This is a good approximation at the local scale (10–100 km) considered in this study, except in regions where topography gradients are important. This will induce loss of precision in our analysis of the stress state in the section on the Pyrenees, but the consistency of the obtained results makes us believe that the inaccuracy is not large enough to invalidate the method, even in this case.

Due to the elastic thickness of the lithosphere, the assumption of isostatic compensation of the lithosphere is not valid at the local scale. Nevertheless, even if the spatial resolution of our study is around 10 km, the spatial extension of most of the regions presenting homogeneous local stress behavior is found to be of larger dimensions, suggesting that the related geoid anomalies have an extension reaching that corresponding to isostatic compensation (100–200 km for a lithosphere with an elastic thickness of 10 km).

STRESSES IN THE LITHOSPHERE

Figure 2A shows, in 1-D, the geoid height (in a relative scale) caused by an arbitrary cylindrical density anomaly centered at x = 0 on a homogeneous lithospheric plate of finite extension from x = –a to x = a. The geoid height N(x) is a continuous function of x that has a maximum value N0 at x = 0 and a value Nmin ∼ 0 at the plate limits x = –a and x = a. Figure 2 also shows the horizontal gradient and second derivative of the geoid height. The function N(x) has a shape similar to the geoid height contribution of any single anomaly of any wavelength that can result from lateral variations of the density in the lithosphere and surface topography. For example, it could represent the geoid height along a section across a mountain, like the Pyrenees. For Earth, the geoid represents the sum of a large number of such contributions at different locations inside the plate. By using the principle of superposition, the analysis that follows can be easily generalized.

In the example of Figure 2A, the plate is in extension with a maximal tensile pressure in the column at x = 0, corresponding to the geoid height maximal value. The tensile pressure decreases progressively to zero, together with the geoid height, when x increases along the positive and negative axis up to x = ±a.

The variation in horizontal stress from one column at x to its neighbor at x + dx is equal to the horizontal body force density forumla (Fig. 2A).

Considering development of the function N(x) in Taylor series around x, we have: 
graphic

As a fast signal power decrease is observed in the Western Europe geoid in function of the degree of the development, we limit the Taylor development of the geoid height around x to the first order.

The force balance equation becomes: 
graphic
Equation 11 expresses the increment of stress (positive or negative) in the column of lithosphere at x from the one in the previous column. The first term of the right-hand side represents this stress increment as it would be if there were no variation of the body force density in the column of lithosphere at x, or in other words, if the geoid gradient is locally constant. The second of these two terms represents the local stress increase or decrease caused by the variations of the body force density in the column at x. Hence, a local extensional (compressive) stress is generated in one column of lithosphere when the difference in the stress variation increases (decreases), corresponding to a negative (positive) second derivative of the geoid. From Equation 11, the strength of this stress is given by 
graphic
Note that these internal stresses created in each column of lithosphere are a part of the total stress state in each column of lithosphere as calculated from the GPE or N (Equation 7). They do not contribute to a global extension or compression of the plate (Fig. 2B), because if the plate is elastically homogeneous, 
graphic

When the geoid gradient is not constant, the local stresses are always present independent of the limit conditions imposed on the plate. In our example (Figs. 2A and 2B), the second derivative of the geoid indicates internal local extension between x = –b and x = b, with a maximal extension at x = 0, and local compression outside this central core of the plate. If the ends of the plate are fixed, there is no global extension of the plate, and the local stress is predominant.

In the general case of the three-dimensional (3-D) Earth, a variation of the body force density is created in a column of lithosphere by the local lateral differences of density, as well as of surface topography with the surrounding columns. Mathematically, it is a two-dimensional (2-D) problem that considers the space variations of the two horizontal components of the body force, which are represented by four space derivatives, forming a 2 × 2 tensor Aij: 
graphic
where forumla
The diagonal elements of the matrix represent the force variation along the horizontal axes. A positive value corresponds to an increase of the force along the axis, which means an extensional force; conversely, a negative value means a compressive force. Hence, the sign the divergence of forumla expresses if the column of lithosphere at (x, y) generates a local compression or extension is given by: 
graphic

If the nondiagonal elements of the matrix Aij are zero, this means an absence of shear forces on the planes normal to the x and y directions, and the principal horizontal stresses follow the axis directions. If the nondiagonal elements differ from zero, the principal horizontal stresses are computed as the eigenvectors of the A matrix (Moler and Stewart, 1973).

Depending on the geometry and the amplitude of the structural anomalies at the origin of the spatial variations of the body force density, the geoid gradient can change dramatically in amplitude and direction over short distances, meaning a strong change of the generated stress in neighboring columns of lithosphere, which departs from the simple 1-D case developed in the beginning of this section.

Knowing the values of the two horizontal principal stresses, the value of forumla is calculated using Equation 3. This allows the three principal stresses to be ordered as σ1 (maximum stress), σ2 and σ3 (minimum stress), and thus the associated tectonic style to be determined and the maximum shear stress value to be calculated as ½(σ1 – σ3).

In our method, we determine the principal horizontal stresses integrated over the whole lithosphere thickness. To obtain the average stress values, we have to divide the integrated stresses by the thickness of the lithosphere. According to the study of Europe by Tesauro et al. (2009), the crustal strength provides a relatively large contribution to the lithospheric strength. They suggested that most of Western Europe is characterized by an elastic thickness between 10 and 20 km, such that a force per unit length of 1010 N/m would correspond to a horizontal stress in the range 0.5–1.0 MPa. In this paper, we only consider orders of magnitude. Consequently, we focus on expressions that do not include the thickness of the lithosphere, as it is not a well-defined parameter everywhere.

In our computations, we consider a surface dimension of the order of 5 km × 5 km for the lithospheric columns, similar to the best geoid resolution. We use the geoid solution EGM2008 (Pavlis et al., 2008), available from the U.S. National Geospatial Intelligence Agency (NGI) with a 2.5 × 2.5 min resolution (http://earth-info.nga.mil/GandG/wgs84/gravitymod/egm2008/). This solution is an optimal combination of space gravity data (Gravity Recovery And Climate Experiment [GRACE] Gravity field and steady-state Ocean Circulation Explorer [GOCE]) and ground gravity anomalies. It is released as a global grid of 8640 × 4321 data points, based on a set of spherical harmonic coefficients up to degree 2190 and order 2159. The coefficients allowed a valid representation down to ∼9 km, but the resolution of the provided 2.5 × 2.5 min grid is twice that value. It is difficult to assess the precision of the geoid data because they are highly dependent on the location. As our application concerns Europe, one of the most monitored areas of the world, we expect that the quality is at the best level of the model, which is at the order of a few centimeters of geoid height, according to validation reports from different authors in Newtons’s Bulletin (2009).

STRESSES AND TECTONIC REGIMES AT THE PASSIVE MARGINS IN THE BAY OF BISCAY

At passive margins, the lower density of continental crust relative to that of the oceanic crust causes tensional stress in the continental part of the margin and compression stress in the oceanic lithosphere (Bott, 1971; Artyushkov, 1973). We study here the passive margin in the Bay of Biscay (Figs. 1, 3, and 4). On the map representing the divergence of the gradient of the geoid anomaly (Fig. 3), the passive margin is associated with a zone with sources of positive divergence on the continental side corresponding to lithosphere under extension, and of negative divergence on the oceanic side associated with lithosphere under compression, as expected from Stein et al. (1989). The orientation of the maximum horizontal principal stress σH is perpendicular to the margin for the sources in the oceanic part near the margin, confirming the expected compression perpendicular to it (Fig. 4). It is parallel to the margin in the continental part, confirming the extension perpendicular to the margin in this side of the margin.

This example, considering a relatively simple and well-understood structural configuration, illustrates the validity of our approach for the interpretation of stress pattern sources caused by lateral heterogeneities and topography in plate interior lithosphere. For the Bay of Biscay margin, our analysis shows that the dimension of the strips of compression and dilatation near the margin is of the order of 50 km. We also found that the shear stress source is the most significant in the lithosphere columns near the top and the bottom of the continental slope, where the value of the integrated shear stress is of the order of 1.2 × 1011 N/m. As the thickness of the elastic lithosphere is 50 km in the oceanic part and 25 km in the continental part of the margin (Tesauro et al., 2009), this corresponds to a shear stress of 2.4 MPa at the bottom and 4.8 MPa at the top of the continental slope. Note that the geoid resolution is sufficient to retrieve the short-wavelength signal despite the fact that this signal is nearly perpendicular to the high-amplitude large-scale geoid variations (Fig. 1).

COMPARISON OF STRESSES FROM THE SECOND SPATIAL DERIVATIVE OF THE GEOID AND THE WORLD STRESS MAP IN WESTERN EUROPE

We apply our method to the case of Western Europe. The divergence of the horizontal force in the lithosphere and the direction of σH calculated from the second spatial derivatives of the geoid (Fig. 5) indicate that the local stress pattern is dominated by relatively short-wavelength variations (a few tens to a few hundreds of kilometers) of the stress regime, with adjacent regions under compression (in blue color) and extension (in red color), which can be superimposed onto the complex geological structures of the region. This is illustrated by the specific case of the Pyrenees in the next section.

In this section, we compare the directions of σH obtained from the second spatial derivatives of the geoid with the ones from the World Stress Map (http://www.world-stress-map.org) (Fig. 6A). Since our stress in each lithosphere column is determined by computing second derivatives from geoid heights at the eight neighboring columns, we compare our results with the WSM by considering the stress state in the column corresponding to the WSM’s data as well as in the eight neighboring columns. In our analysis, we used WSM data on the area extending from 10°W to 10°E in longitude and from 42°N to 53°N in latitude. More than 80% of the 958 available data points are originated from earthquake focal mechanisms (542) and boreholes (254). This data set includes few A-quality and B-quality data (σH measured to within ±15° and ±20°), and near half of the data are D-quality (σH measured to within ±40°). E-quality types were rejected (σH could not be determined or had an uncertainty greater than ±40°).

We note a fair agreement between the two data sets (Fig. 6A). To substantiate this, we performed a statistical test using the null hypothesis that there is no relation between our σH direction estimates and the ones from the WSM, and that the observed agreement is purely coincidental. Under that hypothesis, we randomly generated 100,000 samples of same size as our results based on the second spatial derivative of the geoid at the 958 locations in Western Europe where WSM data exist, each sample being composed of 958 sets Ãj of nine orientations.

Since the nine values from each set from our results have some physical link with each other, comparing them with nine random values distributed uniformly between 0° and 360° would not be a fair comparison. Thus, considering that the distribution of the nine values of each set Ãj is consistent with a normal distribution of standard deviation σi, we drew each set Ãj with a normal distribution of expected value forumla and standard deviation forumla. Each forumla value was randomly taken from one of our 958 observed σi; this allowed the standard deviation forumla to be consistent with our observed σi. On the other hand, each expected value µj was drawn from the WSM values, which is a stronger test than taking one of our mean values.

For each sample and set, we took the minimum angular distance between the nine randomly generated orientations and the WSM ones (modulo 180°) of the nine values, and made the sum of the square of this distance over the 958 sets. At the end of this computation, we created a probability function from the 100,000 values of sums of squared misfits for comparison with our results. We obtained 1.3591e+006 from our σH directions determined from the second spatial derivatives of the geoid, which is smaller than the samples obtained randomly, ranging between 1.4882e+006 and 2.0049e+006 (their distribution is shown on Fig. 7). Consequently, we conclude that our estimate of the agreement in σH direction between our results and the WSM is certainly not the result of random fluctuation of space variable quantities.

Figure 6B presents the frequency distribution of the local maximum horizontal principal direction σH inferred from the second spatial derivative of the geoid at the locations where the WSM provides information. Two main orientations, at 45° (NE-SW) and 135° (NW-SE), are evidenced. These two σH directions interchange in parallel with variations of the tectonic style. They correspond to the two horizontal principal stress directions acting at larger wavelength (200 km or more) in Western Europe (Heidbach et al., 2007, 2010), with the maximum (σH) being oriented NW-SE and the minimum (σh) NE-SW.

The NW-SE direction corresponds to the average direction of the geoid gradient in Western Europe, and hence to the direction of the force acting at large wavelength (several thousands of kilometers) in the lithosphere. This gradient results from the location of the geoid maxima at the Mid-Atlantic Ridge on the western and northwestern side from the considered plate interior and the Alps-Pyrenees system in the S-SE direction (Fig. 1). To explain the bimodal σH direction distribution of the local stress tensor inferred from the second spatial derivative of the geoid, consider 1-D changes of the geoid gradient. It can be generalized in a similar way to the full horizontal case. As the long-wavelength geoid gradient remains oriented NW-SE while the local gradient changes, a local increase of the gradient, corresponding to local sources of extensional force in the NW-SE direction, must be associated with a local decrease of the gradient, corresponding to local sources of compressive force in the same direction, and vice versa. Hence, when the geoid gradient increases, the NW-SE direction reflects the main direction of local extension, and σH is oriented NE-SW. When the geoid gradient decreases, the NW-SE direction reflects the local main direction of compression σH.

The WSM data set (Fig. 6C) includes some data with a σH oriented NE-SW, but this orientation is far less evident than in the local stresses inferred from the second spatial derivative of the geoid (Fig. 6B); 16.7% of the data present a NE-SW σH direction, ranging between 15° and 75° (45° ± 30°), to be compared to the 55.5% of σH oriented around the NW-SE direction (see introduction). To complete the information of Figure 6A, we analyzed the misfit in function of the σH direction from the WSM by ranges of 30°. When σH from the WSM is oriented around the NE-SW direction, it agrees relatively well with the one from our computations, which means that at these locations, the local stress source inferred from the geoid appears as the strongest contribution to the stress state. When σH from the WSM is oriented around NW-SE, the misfit suggests that for a part of the locations where our local stress source is NE-SW, the actual stress state has σH oriented NW-SE. In Western Europe, this direction is parallel to the relative plate motion of the Africa plate with respect to Eurasia plate, and thus at locations where the local stress source is smaller than the far-field stress caused by the collision of the two plates, σH in the NW-SE direction is enhanced. Another possible, but more speculative, contribution to the predominance of the NW-SE σH direction in the WSM data set relative to the bimodal directions inferred from spatial variations in the geoid is anisotropic or variable strength of the continental lithosphere.

The σH direction only provides a part of the stress state characteristics, and using the σH directions from the WSM does not allow a full comparison with our results. The WSM data set includes 542 focal mechanisms; 172 are normal faulting, 40 are normal strike-slip faulting, 72 are thrust faulting, 18 are thrust strike-slip faulting, and 240 are strike-slip faulting events. We compared the tectonic style from our results and the focal mechanism from the WSM. We found that they agree in 71% of the cases; a null hypothesis test showed that such a high value could not be obtained just by chance (the simulation gives a 95% confidence interval of [40%,49%] for the score under the null hypothesis). Considering that the far-field stress field component is supposed to be compressive, even if we suppose that all the strike-slip, thrust strike-slip faulting, and thrust faulting agree with this stress, this corresponds to only 60% of the mechanisms. Therefore, the local stress sources better explain this information than the expected far-field stresses.

To completely address the question, it is necessary to obtain directly a comparison of the different characteristics of the stresses and not only the σH direction and the faulting style. When fault-plane solutions are available, it is possible to define a misfit function as the difference between the slip directions in each of the nodal planes with the shear stress acting in these planes due to the deviatoric stress computed from the second spatial derivative of the geoid. This will be developed in the next section in the case of the Pyrenees, which is an interesting region because a simple seismotectonic model seems inadequate to explain the complexity of earthquake fault-plane solutions (Delouis et al., 1993; Goula et al., 1999; Rigo et al., 2005).

STATE OF STRESS AND EARTHQUAKE FAULT-PLANE SOLUTIONS IN THE PYRENEES

Figure 8 shows the tectonic style and the σH direction deduced from the geoid in the Pyrenees. Focal mechanisms of 44 earthquakes published by Delouis et al. (1993), Goula et al. (1999), Souriau and Pauchet (2001), Rigo et al. (2005), and Chevrot et al. (2011) (Table 1) are also shown. Nine σH directions are represented by small boxes on each earthquake mechanism. They correspond to the directions computed for the column of lithosphere in which the earthquake hypocenter is located and the eight adjacent columns. Unfortunately, the earthquake repartition is not spatially homogeneous, but the spatial variations of the earthquake fault-plane solutions indicate short-wavelength spatial changes of the stress pattern, which is also one of the main characteristics of the stress calculated using the second spatial derivatives of the geoid.

From the tectonic style deduced from the second spatial derivatives of the geoid (Fig. 8), three different quasi–east-west–oriented zones are found. North of the North Pyrenean Frontal thrust, the stress source field is mainly compressive, with a σH direction more or less perpendicular to the North Pyrenean Frontal thrust and a predominance of reverse faulting tectonic style. South of the North Pyrenean Frontal thrust, an area centered on the France-Spain border and corresponding to the North Pyrenean zone and the axial Pyrenean zone, it is mostly characterized by an extensional stress field with an average σH direction parallel to the North Pyrenean Frontal thrust, but the map also indicates small regions characterized by a reverse tectonic style. Most of the seismic activity of the Pyrenees occurs in this zone. The third zone includes the southern part of the map and crosses the South Pyrenean thrust. It appears more complex than the two other ones, with more pronounced spatial variation in the tectonic styles and in the σH direction.

As our stress pattern in each lithosphere column is determined by computing derivatives from geoid heights at the eight neighboring columns, the coherence analysis of fault-plane solutions in the Pyrenees with our stress evaluations is done by considering the stress state in the column corresponding to the earthquake focus but also that in these eight neighboring columns. This corresponds to an uncertainty in the stress pattern location averaging 15 km. For each focal mechanism, we compute the angle (angle 2 in Table 1) between the slip direction of the two nodal planes and the resolved shear stress direction from our deviatoric stress tensor in each of these planes. As indicator of the misfit between the fault-plane solutions and our stress field, we considered the smallest of the 18 determined values (two per lithospheric column). We considered this indicator as a better misfit function than the comparison of the σH direction of our local stress sources with the one obtained using the Zoback (1992) criteria (angle 1), because angle 2 considers the full information of the stress tensor. Figure 9 shows the misfit distribution for the 44 earthquakes considered. Most of the focal mechanisms agree with at least one of the nine stress states from the geoid associated with the earthquake location. The average value of the misfit is 2°, and the standard deviation is 45°. There are earthquakes located in areas where the stress state from the geoid and the associated σH direction change abruptly over a distance comparable to our smallest investigated wavelength of 15 km. This is particularly true for many of the smallest-magnitude earthquakes in our data set. Hence, these data are useless to assess whether our stress field explains the earthquake mechanism, because their small numbers do not allow any relevant statistical treatment. Nevertheless, the concordance with at least one at the nine stress state means that our hypothesis cannot be rejected from these data.

Given that small-magnitude earthquakes are not necessarily representative of the stress field acting at the scale of the lithosphere, we analyze the data from the three recent strongest earthquakes (1, 2, and 8 in Table 1; Fig. 8), with a magnitude greater than 5.0. These earthquakes are likely to represent the deformation field at the regional scale. For them, the nine σH orientations are comparable. The two published mechanisms for the 1967 Arette event (1) are consistent with the WNW-ESE σH direction. The location of the event in a small region where the predicted tectonic style is compressive agrees with the reverse mechanism solution. The strike-slip solution is equally possible, but more discordant with our predicted stresses. The 1980 Arudy (2) as well as the 1996 Saint-Paul de Fenouillet (8) earthquake fault-plane solutions are in good agreement with the tectonic style and σH direction predicted by the local stress sources. Note that only considering σH direction, earthquakes 1 and 2 are consistent with the regional stress pattern, but the difference in their faulting style is well explained by the local stress source. There is also a similar agreement for four (3, 4, 6, 7) of the six earthquakes with 4.4 ≤ M ≤ 5.0. Hence, considering the fault-plane solution of these nine earthquakes and the associated σH direction, seven are explained by the stress resulting from the second spatial derivatives of the geoid, whereas five are explained by the NW-SE far-field stress field. This result suggests that the stress pattern from the local sources better explains earthquake mechanisms in the Pyrenees than the NW-SE far-field stresses.

This conclusion is also supported by the occurrence of most of the natural earthquakes in the columns of lithosphere where the generated maximum shear stress is high (Fig. 10). This is particularly evident in the central part of the Pyrenees between the South Pyrenean thrust and the North Pyrenean Frontal thrust. A parallel can be established with the 2-D case of Figure 2A, roughly corresponding to a geoid transect perpendicular to the Pyrenean Mountains. The region in extension is located between x = ±b, which is equivalent to the central part of the Pyrenees, and where the generated extensional stress is higher than the generated compressive stress outside this limit, as indicated by the curve of d2Ndx2. The maximal depth-integrated local shear stress in this central part of the Pyrenees is around 2 × 1011 N/m. By comparison, the maximal GPE difference between the Pyrenees and its foreland is around 1012 N/m, corresponding to a 4 m difference of the geoid height (Fig. 1).

North of the North Pyrenean Frontal thrust, the local stress is compressive (Figs. 5 and 8), with a maximum shear stress weaker in a large part of the area than in the central zone of the Pyrenees (Fig. 10). There are fewer earthquakes in this region, and most of them occur just north of the North Pyrenean Frontal thrust or more to the west, where the maximum shear stress appears to be the strongest for the area. Few earthquakes occur in the white part of the map, corresponding to small maximal shear stress (Fig. 10). Induced seismic activity is generated in the Lacq deep gas field (Fig. 10); to explain this activity, Segall et al. (1994) suggested a compressive regional stress, which corresponds to the local stress resulting from our computations. Unfortunately, there are no mechanisms available for natural earthquakes in this region to confirm this hypothesis and establish the stress characteristics.

CONCLUSION

We developed, tested, and applied to Western Europe a method to evaluate locally the depth-integrated stress sources within the lithosphere. These local stresses represent the additional body force density created in a column of lithosphere by the lateral differences of density, as well as of surface topography with the surrounding columns. Our method implies the assumption of isostatic compensation of the lithosphere. Due to the elastic thickness of the lithosphere, it is not valid at the spatial resolution of our study, which is around 10 km. Nevertheless, the spatial extension of most of the regions presenting homogeneous local stress behavior, identified on the map of the body force density divergence, is found to be of larger dimensions (Fig. 5). In this case, the stress contribution generated in our small columns of lithosphere represents a part of a larger-wavelength stress associated with geological lateral heterogeneities or (and) topographies that are of sufficient dimension to be isostatically compensated. Thus, they are well representative of this stress generated regionally in the lithosphere. However, in regions where surface relief is pronounced or (and) where short-wavelength density heterogeneities are present in the lithosphere, our results identify very small regions of a few lithospheric columns with a stress different from that in the surrounding columns and should be considered with caution. A solution to overcome this problem would consist of filtering the geoid to get rid of small-scale surface relief or density heterogeneity features. In this study, we prefer to keep our original small-scale resolution, first, because it allows us to better delineate the limits of the regions where extension or compression are generated and compare them with the geological structures, and second, because our interest is in evidence for stress heterogeneities on faults with maximal dimensions of 10–30 km, the limit of the spatial resolution of our method and the potential sources of magnitude 5.5–6.5 earthquakes, the largest known events in Western Europe. It is also worth noting that our method intends to evaluate stresses integrated through the whole thickness of the lithosphere, which corresponds to a low-pass spatial filtering that already suppresses the shorter wavelengths.

In Western Europe, our results indicate that the stress pattern from the local sources is dominated by short-wavelength changes (few tens to hundreds of kilometers) in the tectonic style and the σH direction. It appears to be directly related to the regional geological structures, which is not really a surprise given that these structures are at the origin of the mass heterogeneities in the lithosphere.

The comparison of the σH orientations and tectonic style from the local sources with the ones from the World Stress Map (WSM) data set indicates that the local stress sources can be representative of the deviatoric stress state in some regions, but not everywhere. This is because the σH direction NW-SE is enhanced in the WSM by comparison to the NE-SW direction, which is as well distributed as the NW-SE direction in the local sources. Because the NW-SE σH direction is well represented in our results and the WSM data, it is not easy to identify the regions where the local stress sources are preponderant. Nevertheless, it is important to note that the local stress sources explain 71% of the faulting styles for the earthquake fault-plane solutions in the WSM. This is better than the 60% agreement given by the classical compressive NW-SE far-field stress model resulting from the collision between Africa and Europe, and from the opening of the Atlantic Ocean. Heidbach et al. (2007, 2010) already pinpointed that local stresses can control the stress pattern at different spatial scales in Western Europe. Our study not only confirms their result, but it also provides a method with which to evaluate these local stresses, whatever their wavelength.

In the Pyrenees, we obtained an agreement for 44 earthquakes between their fault-slip directions and the direction of shear stress from the local sources acting on the associated fault planes. Unfortunately, the number of data in the data set and its poor spatial repartition are not sufficient to prove by a statistical criterion that local stress sources explain the earthquake mechanisms. However, considering the fault-plane solution and the associated σH direction of the nine earthquakes in our data set with M greater than 4.5, seven are explained by these local stress sources, whereas five are explained by the far-field stresses. This result suggests that the stress pattern from the local sources could better explain earthquake mechanisms in the Pyrenees than the classical model. Our stress field from the second spatial derivative of the geoid is compatible with the extensional stress field evidenced by Chevrot et al. (2011) in the central Pyrenees.

These comparisons of the local stress sources from the second spatial derivative of the geoid with the WSM data set and stress characteristics from earthquake fault-plane solutions, which are supposed to represent the stress state in the lithosphere, are motivation for further using the geoid in this manner. As the geoid is a well-determined surface everywhere on the Earth, our method allows the stress source pattern to be evaluated everywhere in continental and oceanic plate interiors worldwide. The method does not require modeling of the lithosphere mechanical parameters, which is necessary when modeling the GPE. Future studies should be devoted to constraining the order of magnitude for which the use of the geoid instead of the GPE is valid by comparing our results inferred from the second spatial derivatives of the geoid with the ones using the GPE obtained by modeling the lithosphere elastic properties.

In 2008, the European Space Agency GOCE mission was launched: With nearly 2 yr of data, it has already provided a global geoid undulation model of unprecedented precision and resolution (200 km) (Pail et al., 2010). Our study suggests an application of those data to investigate the stress in the lithosphere. On the other hand, the small-scale signal that we identify shows the complementary nature of high-precision global space data with high-resolution terrestrial gravity data (Panet et al., 2011).

We thank Alexis Rigo for the fruitful discussion on the earthquake activity in the Pyrenees and for providing us the paper of Chevrot et al. (2011) before its publication. D. Kusters is supported by FRIA (Fund for the Research in the Industry and Agriculture) contract FC 86440. The work of O. de Viron and M. Van Camp constitutes Institut de Physique du Globe de Paris contribution number 3351. We thank the editor, Raymond M. Russo, the reviewers, James Conder and Laurent Husson, and a third anonymous reviewer for their valuable comments, which helped us to improve the original manuscript. This study benefited from the support of the Campus Spatial, and from CNES (Centre national d’études spatiales) through the TOSCA program, as an exploitation of the GRACE mission data. The U.S. National Geospatial Intelligence Agency (NGI) is acknowledged for making the geoid data available.