Lithospheric and sublithospheric deformation under the Borborema Province of northeastern Brazil from receiver function harmonic stripping

The depth-dependent anisotropic structure of the lithosphere under the Borborema Province in northeast Brazil has been investigated via harmonic stripping of receiver functions developed at 39 stations in the region. This method retrieves the first (k=1) and second (k=2) degree harmonics of a receiver function dataset, which characterize seismic anisotropy beneath a seismic station. Anisotropic fabrics are in turn directly related to the deformation of the lithosphere from past and current tectonic processes. Our results reveal the presence of anisotropy within the crust and the lithospheric mantle throughout the entire province. Most stations in the continental interior report consistent anisotropic orientations in the crust and lithospheric mantle, suggesting a dominant northeast–southwest pervasive deformation along lithospheric-scale shear zones developed during the Brasiliano–Pan-African orogeny. Several stations aligned along a northeast–southwest trend located above the (now aborted) Mesozoic Cariri–Potiguar rift display large uncertainties for the fast-axis direction. This non-azimuthal anisotropy may be related to a complex anisotropic fabric resulting from a combination of deformation along the ancient collision between Precambrian blocks, Mesozoic extension and thermomechanical erosion dragging by sublithospheric flow. Finally, several stations along the Atlantic coast reveal depth-dependent anisotropic orientations roughly (sub)perpendicular to the margin. These results suggest a more recent overprint, probably related to the presence of frozen anisotropy in the lithosphere due to stretching and rifting during the opening of the South Atlantic.

1 Introduction

Understanding intraplate deformation and its relationship to deep geodynamic processes such as sublithospheric flow is critical for improving our understanding of the evolution of continents. The Borborema Province in northeastern Brazil, for instance, has witnessed several cycles of deformation, as well as recurrent episodes of intraplate volcanism and uplift, during its geological history. Brasiliano–Pan-African deformation is well represented through the network of shear zones that pervade the Borborema Province (Vauchez et al., 1995; Neves et al., 2000). These shear zones separate several tectonic terrains of the Paleoproterozoic Era and Archean Eon that amalgamated and/or were reworked during the orogeny (Jardim de Sá et al., 1992; Cordani et al., 2003). Major Neoproterozoic shear zones thus constitute inherited structures that could have influenced the geometry of subsequent tectonic processes, such as the opening of the South Atlantic Ocean (Tommasi and Vauchez, 2001; Kirkpatrick et al., 2013). Also, the current topography of the Borborema Plateau and the Sertaneja Depression may have resulted from a combination of ongoing deep processes, such as edge-driven convection in the asthenospheric mantle and/or stretching and thinning of the lithosphere during Mesozoic times (de Oliveira and Medeiros, 2012; Almeida et al., 2015).

Recent seismological studies from receiver functions (Almeida et al., 2015; Pinheiro and Julia, 2014; Luz et al., 2015a; Luz et al., 2015b), ambient noise (Dias et al., 2014) or P-wave tomography (Simões Neto et al., 2019), and SKS splitting (Bastow et al., 2015) have greatly contributed to further our understanding of the relationships between inherited Precambrian structures, Mesozoic extensional processes, and episodes of post-breakup volcanism and uplift. However, several tectonic and geodynamic questions remain unanswered. In particular, the presence of the Meso-Cenozoic Macau–Queimadas volcanism (MQA, Fig. 1) – which does not present a clear age progression – remains unclear. Moreover, SKS waves showed little to no evidence of splitting in the continental interior (Bastow et al., 2015), which is difficult to comprehend given the complex tectonic and deformational history of the province.

Here, we determine depth-dependent anisotropy in the Borborema lithosphere (crust and mantle) from harmonic analysis of receiver functions. Our results confirm that SKS splitting at coastal stations is dominated by fossil anisotropic fabrics in the lithospheric mantle, likely originating from Mesozoic extension. In the continental interior, receiver function stripping reveal fast-axis orientations consistent with major regional shear zones, suggesting their continuation at depth into the lithospheric mantle. Our results also show the presence of non-azimuthal anisotropy above stations located along the now aborted Cariri–Potiguar rift. This non-azimuthal anisotropy is likely related to complex fossil anisotropic fabrics, resulting from a combination of deformation along the ancient collision between Precambrian blocks, the Mesozoic extension, and thermomechanical erosion/mantle dragging by sublithospheric flow.

2 Geological setting

Figure 1 The main geological features of the Borborema Province superimposed on a topographic map of northeastern Brazil. Black and grey lines correspond to major shear zones (SZ), and red dashed lines correspond to the volcanic alignments of Fernando de Noronha–Mecejana (FNMA) and Macau–Queimadas (MQA). The Borborema Plateau boundaries are indicated in blue.

The Borborema Province formed during the Neoproterozoic Brasiliano–Pan-African orogeny (600–580 Ma), as a result of the collision between the São Luiz–West Africa Craton to the north and the São Francisco–Congo Craton to the south (Jardim de Sá et al., 1992; Cordani et al., 2003). Thus, it represents the west central portion of a larger Neoproterozoic belt that resulted from the assembly of the Gondwana supercontinent.

The basement of the Borborema Province comprises mostly gneisses and migmatitic rocks of Paleoproterozoic age, and small Archean nuclei, overlain by Neoproterozoic metasediments formed during the Brasiliano orogeny (Neves, 2003). This basement is affected by an extensive network of Neoproterozoic shear zones oriented east–west and northeast–southwest (Fig. 1). These shear zones are major structures several hundreds of kilometers long and tens of kilometers wide (Vauchez and da Silva, 1992) that can be traced into the African continent in paleogeographic reconstructions (Arthaud et al., 2008). The Borborema shear zones were activated in high temperature and high-to-low pressure conditions and are associated with a strong production of magmas from both crustal and mantle sources (Vauchez et al., 1995). The shear zone network can be split into two domains: a western domain of northeast-striking strike-slip faults, and an eastern domain of more sinuous, discontinuous east–west-striking shear zones (Vauchez et al., 1995). These two domains could be related to two discrete collisional events with the Parnaíba Block to the west and the São Francisco Craton to the south, respectively, which forced northeast extrusion of the province at the end of the Neoproterozoic (Araújo et al., 2014).

The geodynamic evolution of this basement and the significance of these shear zones is still debated, and two main models have traditionally been proposed. On the one hand, the accretionary model proposes that the Borborema Province is comprised of several Paleoproterozoic small continental fragments that aggregated along the shear zones, which then constitute lithospheric-scale suture zones separating independent tectonic blocks (Cordani et al., 2003; Van Schmus et al., 2011; Araújo et al., 2014). The number of independent terrains is unclear, but there is a general consensus in arranging them into five major Precambrian domains: (i) Médio–Coereau, in the northwestern most tip of the province; (ii) Ceará, between the Sobral–Pedro II and Jaguaribe–Tatajuba shear zones; (iii) Rio Grande do Norte, immediately east of the Ceará domain; (iv) transversal or central, between the Patos and Pernambuco lineaments; and (v) southern, immediately north of the São Francisco Craton (Fig. 1). On the other hand, it has been suggested that the Borborema Province has been a single unit since 2.0 Ga and that the shear zones recorded intracontinental supracrustal deformation during the Brasiliano orogeny (Tommasi et al., 1995; Vauchez et al., 1995; Neves, 2003). In this later model, micro-plate amalgamation largely predates the Brasiliano orogeny, which would have only partially reworked preexisting plate structures in the Borborema Province.

Figure 2 Topographic map of northeast Brazil with locations of broadband and short-period stations considered in this study. Stations are color-coded by network: stations from the RSISNE network are represented by dark blue, INCT-ET by green, Milenio by red, GSN by pink, BLSP by yellow, BODES by light blue and “other” networks by grey (see legend). Only the selected stations have been named and are represented using black contours.

During the Mesozoic, opening of the South Atlantic separated the Borborema Province from its African conjugate (de Matos, 1992). Continental rifting resulted in significant crustal thinning in the region (Santos et al., 2014; Lima Neto et al., 2013; Luz et al., 2015b), forming both marginal (e.g., Ceará, Potiguar, Pernambuco–Paraíba and Sergipe–Alagoas) and intracontinental (e.g., Araripe and Tucano) sedimentary basins (Fig. 1). The post Gondwana breakup evolution of the Borborema Province is characterized by recurrent magmatism (Knesel et al., 2011) and postulated episodes of uplift in the Borborema Plateau (Morais Neto et al., 2009; de Oliveira and Medeiros, 2012) and the Araripe Basin (Assine, 2007; Marques et al., 2014). Intraplate volcanism is characterized as small-volume, long-lived and mainly alkalic in nature (Knesel et al., 2011). It is arranged along two main linear alignments of Mesozoic–Cenozoic volcanic rocks (Fig. 1): the Macau–Queimadas Alignment (MQA), mostly onshore and approximately trending in the north–south direction; and the Fernando de Noronha–Mecejana Alignment (FNMA), mostly offshore and trending in the east–west direction (Mizusaki et al., 2002; Knesel et al., 2011, and references therein). The MQA displays K/Ar ages ranging from 80 to 30 Ma (Mizusaki et al., 2002) and 40Ar/39Ar ages ranging from 93 to 7 Ma (Knesel et al., 2011) without a clear age progression, whereas the FNMA displays progressive K/Ar and 40Ar/39Ar ages from the Fernando de Noronha archipelago (22 to 2 Ma) to the west (Knesel et al., 2011, and references therein), to the Mecejana volcanism (34 to 26 Ma) to the east (Mizusaki et al., 2002).

Cenozoic uplift along the northeastern Brazilian margin was inferred from relative dating of elevated sediments of the Serra dos Martins Formation, absolute dating from apatite fission-track analysis of granitic-gneissic and sedimentary samples and geomorphological studies (Morais Neto et al., 2009; de Oliveira and Medeiros, 2012; da Nóbrega et al., 2005; Nogueira et al., 2015). Although some sort of tectonic uplift and/or inversion of the Araripe Basin seems to be widely accepted (Marques et al., 2014; Peulvast and Bétard, 2015; Garcia et al., 2019), uplift in the Borborema Plateau is more debated. de Oliveira and Medeiros (2012) argue that this tectonic event is linked to Cenozoic mafic underplating and isostatic uplift due to a small-scale convection cell at the edge of the continent, which might have also been responsible for the surface volcanism (Knesel et al., 2011). The hypothesis of a thin layer of mafic underplate seems to be consistent with recent receiver functions observations south of the Patos Lineament (Luz et al., 2015a, b). However, Luz et al. (2015b) debated the time of emplacement of such mafic cumulates. These authors proposed that this mafic layer would be part of the original Proterozoic crust, and that the southern Borborema Plateau should be regarded as a high-standing, rheologically strong block surrounded by stretched and delaminated crust. The stretching model seems to have been confirmed by a rheological contrast along the Patos Lineament postulated from seismic P-wave tomography (Simões Neto et al., 2019). The tomographic study also identifies an asthenospheric low-velocity channel trending northeast–southwest under the center of the province, which is interpreted as resulting from lateral flow from a distant mantle plume. Such asthenospheric flow might represent the source of Meso-Cenozoic intraplate volcanism in northeastern Brazil.

3 Data and methodology

3.1 Seismic data

Seismic data for this study were obtained at 75 seismic stations in northeastern Brazil. These stations belong to a variety of seismic networks, both permanent and temporary. The Rede Sismográfica do Nordeste (RSISNE) consists of 19 broadband stations equipped with RefTek 151-120 sensors feeding RT-130 digitizers (24-bit) sampling at 100 Hz, with an interstation spacing of about 250 km and a network aperture of ∼800 km. The RSISNE network has been in operation since 2011 and was initially funded by the national oil company Petrobras. The Instituto Nacional de Ciência e Tecnologia em Estudos Tectônicos (INCT-ET) network, consists of 7 broadband stations and 22 short-period stations. The broadband stations were arranged along an approximately 1000 km long line with interstation spacing of about 100 km. They were equipped with STS-2.5 Streckeisen sensors and Q330 data loggers (24-bit) sampling at 100 Hz. The 22 short-period stations were equipped with L4A-3D Sercel sensors (2 Hz cutoff frequency) and 24-bit RT-130 digitizers sampling at 100 Hz. They were in operations between 2011 and 2012 and recorded continuously for between 6 months and 1.5 years. The INCT-ET network was funded by the Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq). Up to six broadband stations operated during the period from 2007 to 2009 under the “Institutos do Milênio” project. These stations were equipped with KS-2000 Geotech sensors and Geotech digitizers sampling continuously at 100 Hz. They recorded continuously for periods ranging from 6 months to 2 years, with two of them still currently in operation. This network was also funded by CNPq. Station RCBR belongs to the Global Seismographic Network (GSN). This station has been recording since March 1999 using a Guralp CMG-3T sensor (replaced in July 2004 by a STS-2 Streckeisen sensor), which always feeds a Q330 data logger and samples continuously at 40 Hz. Seven broadband stations belonging to the broader Brazilian Lithosphere Seismic Project (BLSP) (Assumpção et al., 2004) operated for 1.5 to 3.0 years in northeastern Brazil. They were equipped with either Guralp CMG-3T or STS-2 Streckeisen sensors and 24-bit RT-130 digitizers sampling at 100 Hz. Finally, 11 broadband stations deployed during the Borborema Deep Electromagnetic and Seismic (BODES) experiment were installed along an approximately north–south line crossing the Araripe Basin. They were equipped with RefTek 151-120 sensors and RT-130 digitizers. The stations were in operation between 2015 and 2017 and recorded continuously for ∼2 years at 100 Hz. Further details regarding the 75 seismic stations are given in the Table S1, and their geographical location is displayed in Fig. 2.

Figure 3 Location of earthquakes (blue circles) used for receiver function analyses, occurring at epicentral distances between 30 and 90^∘ (red lines) and with a magnitude of M_b≥5.0. The yellow triangle corresponds to the seismic network presented in Fig. 2.

3.2 Receiver function processing and migration

Receiver functions were computed for the 75 stations comprising the combined network for northeast Brazil. Most of the receiver function estimates were developed by Luz et al. (2015a) in order to investigate lateral variations in crustal thickness and the bulk Vp ∕ Vs ratio across the Borborema Province. This dataset was later utilized by Almeida et al. (2015) to study the crustal architecture of the region from receiver function common conversion point (CCP) stacks. We developed 1400 new receiver function estimates at 11 temporary stations from the BODES network, following the same procedure as Luz et al. (2015a). The receiver function approach aims at retrieving P-to-S converted phases within the coda of teleseismic P waves that result from the interaction of the teleseismic P-wave-front with crustal and upper mantle discontinuities under the recording station (Langston, 1979, 1977) in order to produce estimates of the depth of the discontinuities. The converted phases are detected by deconvolving the radial and tangential components of the teleseismic waveforms by the corresponding vertical component (Ammon, 1991). This operation removes the effects of the source time function, near source propagation and instrumental response from the seismograms, leaving the signature of propagation local to the receiver.

Figure 4 Example of stacked receiver functions represented by back-azimuth bins of 10^∘ at station PFBR. Grey numbers correspond to the number of stacked receiver functions.

The main processing steps involved in the development of the receiver function estimates are summarized below, and further details regarding computational and quality control procedures can be found in Luz et al. (2015a). First, we selected seismic sources with magnitude greater than 5.0 Mb and occurring at epicentral distances between 30 and 90^∘ from the selected stations (see Fig. 3). The corresponding waveforms were then windowed 10 s before and 110 s after the P-wave arrival time, demeaned, detrended, tapered with a 5 % cosine taper, and high-pass filtered above 0.05 Hz to remove low-frequency noise. All waveforms were resampled to 20 Hz, after low-pass filtering below 8 Hz to avoid aliasing. Before deconvolution, the waveforms were additionally low-pass filtered below 1.25 Hz with an acausal Gaussian filter (Gaussian width 2.5). The deconvolution procedure of the vertical component from the radial and transverse components was implemented using the iterative, time-domain procedure of Ligorria and Ammon (1999), with 500 iterations. The deconvolved time series were again filtered with the same Gaussian filter (width 2.5). Percent recoveries of the observed radial component under 85 % were automatically rejected, and the remaining receiver functions were visually inspected for each station to identify and remove outliers.

Prior to implementing the anisotropy analysis, each radial and tangential receiver function was migrated to depth after P to S ray-tracing through the global velocity model ak135-f (Kennett et al., 1995; Montagner and Kennett, 1996). The purpose of the migration is to correct the phase move-out introduced by varying incidence angles among the incoming teleseismic P-wave-fronts, effectively equalizing the receiver function waveforms in the depth domain (Dueker and Sheehan, 1997). Migration before harmonic stripping at individual stations has previously been utilized by Audet (2015), Cossette et al. (2016) and Tarayoun et al. (2017). Similarly, Bianchi et al. (2010), Piana Agostinetti et al. (2011) and Piana Agostinetti and Miller (2014) applied harmonic decomposition on depth-migrated cross sections obtained through CCP stacking of receiver functions. Next, the migrated radial and transverse receiver functions for each station were grouped by back azimuth in 36 nonoverlapping, 10^∘ wide bins, and averaged within each bin. A given station was then selected if it presented two averaged receiver functions (one radial and one tangential) in at least nine bins. This selection criterion ensured a sampling of at least 90^∘ in back azimuth, either continuously or discontinuously, around the station. A back-azimuthal coverage from at least nine bins (each 10^∘ wide) allows the mapping of either half the period for a two-lobed pattern (anisotropy with plunging fast axis of symmetry) or a full period for a four-lobed pattern (anisotropy with horizontal fast axis of symmetry). A total of 39 stations were thus selected for anisotropy analysis. An example of stacked and migrated receiver functions is displayed in Fig. 4.

Figure 5 Example of results obtained at station PFBR using the harmonic stripping method of Bianchi et al. (2010). (a) The harmonic functions obtained by solving Eq. (1) are represented from top to bottom: A(z), first harmonic degree (k=0); B(z), cos (ϕ) term of the second degree harmonic (k=1); C(z), sin (ϕ), term of the second degree harmonic (k=1); D(z), cos (2ϕ) term of the third (k=2) harmonics; and E(z), sin (2ϕ) term of the third (k=2) harmonics. (b) The energy is represented for harmonic degrees k=1 and k=2.

3.3 Estimating depth-dependent anisotropy within the lithosphere

In order to map deformation within the lithosphere, we estimate seismic anisotropy from the harmonic decomposition of receiver functions. The harmonic stripping method is described in Shiomi and Park (2008), Bianchi et al. (2010) and Audet (2015). The method assumes that, at every depth, an ensemble of receiver functions can be expressed as a linear combination of cos (kϕ) and sin (kϕ) terms, where k is the harmonic degree or order, and ϕ is the back azimuth. Shiomi and Park (2008) show that, for anisotropic media, radial and tangential receiver functions display a π∕2k shift for both k=1 and k=2 harmonic degrees; the tangential receiver functions can thus be added to the radial component after applying a phase shift of $+\mathit{\pi }/\mathrm{2}k$ and naturally improve the azimuthal coverage around the station. After the harmonic decomposition is performed, up to five coefficient functions, corresponding to the first three harmonics, are obtained ($k=\mathrm{0},\mathrm{1},\mathrm{2}$). The first harmonic (k=0) represents the isotropic variations from flat interfaces in an equivalent isotropic medium; for this harmonic, the signal is only present in the radial component. If anisotropic structures are present at depth, the second and third harmonics (k=1 and k=2) contain energy with periodicity of 2π∕k. For k=1, a two-lobed periodicity of 2π is either related to the presence of a dipping interface or to an anisotropic layer with a plunging symmetry axis (Maupin and Park, 2007). Two coefficient functions express the projection of this harmonic along the north–south and east–west directions, which correspond to the coefficients multiplying the cos (ϕ) and sin (ϕ) terms, respectively. For k=2, a four-lobed periodicity of π is related to the presence of an anisotropic layer with a horizontal symmetry axis (Maupin and Park, 2007). As for the second degree harmonic, two coefficient functions express the projection of this harmonic degree along the north–south and 45^∘ N directions, corresponding to the coefficients multiplying the cos (2ϕ) and sin (2ϕ) terms, respectively. The harmonic decomposition can be expressed in matrix form (Eq. 1) and solved for the five coefficients for the three harmonic degrees ($k=\mathrm{0},\mathrm{1},\mathrm{2}$) through a singular value decomposition. These harmonics are calculated for every depth within a selected depth-window.

The matrix equation that implements the harmonic decomposition is given by

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\left[\begin{array}{c}{R}_{\mathrm{1}}\left(z\right)\\ \mathrm{⋮}\\ R_n\left(z\right)\\ T_{\mathrm{1}}\left(z\right)\\ \mathrm{⋮}\\ T_n\left(z\right)\end{array}\right]=\\ {\displaystyle }& {\displaystyle }\left[\begin{array}{ccccc}\mathrm{1}& \cos \left({\mathit{\varphi }}_{\mathrm{1}}\right)& \sin \left({\mathit{\varphi }}_{\mathrm{1}}\right)& \cos \left(\mathrm{2}{\mathit{\varphi }}_{\mathrm{1}}\right)& \sin \left(\mathrm{2}{\mathit{\varphi }}_{\mathrm{1}}\right)\\ \mathrm{⋮}& \mathrm{⋮}& \mathrm{⋮}& \mathrm{⋮}& \mathrm{⋮}\\ \mathrm{1}& \cos \left({\mathit{\varphi }}_n\right)& \sin \left({\mathit{\varphi }}_n\right)& \cos \left(\mathrm{2}{\mathit{\varphi }}_n\right)& \sin \left(\mathrm{2}{\mathit{\varphi }}_n\right)\\ \mathrm{0}& \cos ({\mathit{\varphi }}_{\mathrm{1}}+\mathit{\pi }/\mathrm{2})& \sin ({\mathit{\varphi }}_{\mathrm{1}}+\mathit{\pi }/\mathrm{2})& \cos \left(\mathrm{2}\right({\mathit{\varphi }}_{\mathrm{1}}+\mathit{\pi }/\mathrm{4}\left)\right)& \sin \left(\mathrm{2}\right({\mathit{\varphi }}_{\mathrm{1}}+\mathit{\pi }/\mathrm{4}\left)\right)\\ \mathrm{⋮}& \mathrm{⋮}& \mathrm{⋮}& \mathrm{⋮}& \mathrm{⋮}\\ \mathrm{0}& \cos ({\mathit{\varphi }}_n+\mathit{\pi }/\mathrm{2})& \sin ({\mathit{\varphi }}_n+\mathit{\pi }/\mathrm{2})& \cos \left(\mathrm{2}\right({\mathit{\varphi }}_n+\mathit{\pi }/\mathrm{4}\left)\right)& \sin \left(\mathrm{2}\right({\mathit{\varphi }}_n+\mathit{\pi }/\mathrm{4}\left)\right)\end{array}\right]\\ \text{(1)}& {\displaystyle }& {\displaystyle }\times \left[\begin{array}{c}A\left(z\right)\\ B\left(z\right)\\ C\left(z\right)\\ D\left(z\right)\\ E\left(z\right)\end{array}\right],\end{array}$$

where ϕ_i is the back azimuth of the i-th R (radial) and T (tangential) receiver function doublet, A(z) represents the first harmonic coefficient (k=0), B(z) and C(z) are the coefficient functions of the second harmonic (k=1), and D(z) and E(z) are the coefficient functions of the third (k=2) harmonic.

After solving the matrix Eq. (1) within a specific depth-window and calculating the five harmonic coefficients, we search for the presence of anisotropy by inspecting the B(z), C(z), D(z) and/or E(z) terms. If at least one of these component displays nonzero amplitudes, we calculate the energy of the second (k=1) and the third (k=2) harmonic degrees as proposed by Licciardi and Piana Agostinetti (2016):

$$\begin{array}{}\text{(2)}& E_{k=\mathrm{1}}=\sum _{z=\mathrm{1}}^n\left(B(z{)}^{\mathrm{2}}+C(z{)}^{\mathrm{2}}\right)\end{array}$$

and

$$\begin{array}{}\text{(3)}& E_{k=\mathrm{2}}=\sum _{z=\mathrm{1}}^n\left(D(z{)}^{\mathrm{2}}+E(z{)}^{\mathrm{2}}\right).\end{array}$$

These energy functions allow for discrimination between dipping interfaces or anisotropy with a plunging axis of symmetry and anisotropy with an horizontal axis of symmetry. If $E_{k=\mathrm{1}}>E_{k=\mathrm{2}}$, the dominant anisotropy is either a dipping interface or an anisotropic layer with a plunging axis of symmetry. In that case, we rotate B(z) and C(z) in discrete back-azimuth increments α (where $\mathit{\alpha }\in [\mathrm{0},\mathrm{2}\mathit{\pi }]$) and search for the value of α that maximizes B(z), and therefore minimizes C(z). This value of α can be directly interpreted as either the trend of the dip, in the case of dipping interface, or as the trend of the fast axis of symmetry in the case of plunging axis of symmetry. If $E_{k=\mathrm{2}}>E_{k=\mathrm{1}}$, the dominant anisotropy is an anisotropic layer with a horizontal axis of symmetry. In that case, we rotate D(z) and E(z) for each angle increment α (where $\mathit{\alpha }\in [\mathrm{0},\mathrm{2}\mathit{\pi }]$) and search for the value of α that maximizes D(z), and therefore minimizes E(z). This value of α can be directly interpreted as the trend of either the fast or the slow axis of symmetry.

Figure 6 (a) Map of symmetry directions (dark lines) obtained for the crust (0–32 km). When one line is plotted at the station, it represents either the trend of the dip in the case of dipping interface or the trend of the fast axis in the case of plunging anisotropy. When two lines are plotted, they refer to the fast axis and to its perpendicular direction in the case of horizontal anisotropy. Light colors represent 2σ uncertainties estimated from the bootstrap (after resampling 200 times). (b) Same for the lithospheric mantle (32–100 km). Station symbols have been color-coded according to the energy level of the dominant harmonic degree.

An example of harmonic decomposition is shown for station PFBR in Fig. 5. In order to estimate uncertainties, we applied a bootstrap statistical approach by randomly resampling our receiver functions with replacement. We performed the analyses with 200 replications at each of the selected stations. From these 200 values, we estimated the standard error (2σ), which corresponds to the uncertainty in the direction of the fast axis of symmetry. A measurement is considered as unreliable, and then rejected, if the estimated uncertainties are larger than 20^∘.

Table 1 Results of anisotropic symmetry directions for several depth ranges: 0–32, 32–100 and 0–100 km. One direction corresponds to either the trend of the dip in the case of dipping interface or to the trend of the fast axis in the case of plunging anisotropy. When two directions are indicated, they refer to the fast axis and to its perpendicular direction (horizontal anisotropy). Uncertainties were estimated from bootstrap quantification (resampling 200 times at every station).

4 Results

Anisotropy parameters were examined for each station at two depth-window ranges: (1) crust (Fig. 6a), which was assumed to be located at a depth of between 0 and 33 km, in agreement with the 32–40 km range estimated by Luz et al. (2015b) under the Borborema Plateau and 30–33 km under the surrounding basins; and (2) lithospheric mantle, which was taken to be at a depth between 33 and 100 km (Fig. 6b). We assume that the layer with the strongest anisotropy will dominate the results in the case of several anisotropic layers. However, it might happen that results reflect the average value from different anisotropic layers, or from different types of anisotropy in the case of similar anisotropic strength. All results are indicated in Table 1.

Figure 7 Comparison between fast axis of symmetry recorded by SKS waves (red lines) and k=2 harmonics (green lines). SKS-splitting results are from Bastow et al. (2011) and Assumpção et al. (2011). Red lines refer to mean fast-axis orientation (line direction) and delay time (line size) beneath the station. When SKS measurements provide only null measurements, we display black lines which are in the direction of the back azimuth of the recorded event.

An inspection of Fig. 6a reveals that the crust of northeast Brazil presents seismic anisotropy, both within the interior of the continent and along the coast. A number of stations, however, display uncertainties larger than 20^∘. The significance of these large uncertainties are discussed in Sect. 5. Unresolved anisotropic directions within the crust are recorded around longitude −40^∘ for stations nbpb, ar02, ar05, ar06 and nbpn, at the border of the Borborema Plateau (stations nbta, pctv, nbli and caub), and within the Sergipe–Alagoas and Pernambuco basins (stations nban and pcal). The majority of stations that sample clear anisotropic directions display a northeast–southwest to east–west trending axis of symmetry, except stations cs6b (trends ∼ north-northwest–south-southeast), km60 and nbma (trend north–south). We mainly measure anisotropy with 2π-periodicity (k=1) related to a dipping interface or anisotropy with a plunging axis of symmetry, but some stations display clear π-periodic horizontal anisotropy (stations ar09, sabr, pcsa, pcsl, pcja, lp06 and nbmo). We note that, in most cases, even though k=1 harmonics display higher energy contents than k=2 harmonics, both pairs of harmonics display energies with comparable strengths (see Supplement, Fig. S1). For example, station PFBR (Fig. 5) shows clear, nonzero energy levels for both k=1 and k=2 harmonics in the crust between 0 and 33 km.

Figure 6b shows that the lithospheric mantle is characterized by seismic anisotropy throughout the entire province, with the exception of a few stations that display large uncertainties (discussed in Sect. 5). Those include stations within the Parnaíba Basin (trsb), around longitude −40^∘ (ar01, ao05 and nbpn), along the southern portion of the Borborema Plateau (nbta and pcse) and along a northeast–southwest axis located northwest of the Borborema Plateau (nbma, pfbr, nbpa and cs6b). Most anisotropic directions trend northeast–southwest to east–west, with the exception of stations km60 and nbma that show north–south trends. As for the crust, we mainly measure anisotropy with 2π-periodicity (k=1) related to a dipping interface or anisotropy with a plunging axis of symmetry, but some stations display clear π-periodic horizontal anisotropy (stations ar02, km60, lp06, pcal, nban, nbmo, nbpb, pctv and rcbr). Note that stations located along the continental margin show anisotropy within the lithospheric mantle with a fast axis of symmetry that is oblique (stations nbmo, nbpv, nbrf, nban and nbit) or perpendicular (stations nbcl and pcal) to the coast.

In the case where stations recorded anisotropic directions at both crustal and mantle levels, most of them show consistent orientations in the two domains. There are a few instances, nonetheless, in which unaligned orientations for the crust and the lithospheric mantle are observed (ar04, ar09, nbcl, nbit, nbps, nbrf and rcbr).

5 Discussion

A complex combination of lattice-preferred orientation (LPO) and shape-preferred orientation (SPO) could be present in the mantle, although LPO is likely to dominate (Nicolas and Christensen, 1987; Silver, 1996; Mainprice et al., 2000). Fractures and cracks or fine layering, could additionally contribute in the crust. For that reason, our interpretations focus dominantly on mantle anisotropy, consistency of anisotropy within the lithosphere (crust and mantle), and regional-scale trends. Moreover, to avoid a bias related to local features, we refrain from interpreting small-scale variations in anisotropy within the crust.

5.1 Pervasive anisotropy with a (sub)horizontal fast axis of symmetry

As described is Sect. 4, we observe a dominance of 2π-periodicity (k=1) anisotropy in the lithospheric mantle in northeast Brazil, which represents either a dipping interface or anisotropy with a plunging axis of symmetry. However, a close inspection of the energy of the k=2 harmonics also suggests an important contribution from anisotropy with a horizontal axis of symmetry. We were able to replicate this pattern with synthetic receiver functions by assuming anisotropy with a slightly (10 to 15^∘) dipping axis of symmetry (see Supplement, Fig. S2). Note that, in that case, both k=1 and k=2 harmonics display consistent orientations. This is, the orientation inferred from the k=1 harmonic degree is always parallel to one of the two orientations inferred from the k=2 harmonic degree.

Moreover, we noticed that the anisotropic fast axes of symmetry throughout the crust are consistent with those throughout the lithosphere for most of the stations, suggesting a prolongation of crustal structures within the lithospheric mantle. Within the continental interior, the anisotropic orientations are parallel or subparallel to the main east–west to northeast–southwest shear zone directions (stations ar02, nbli or lp06, for example). The consistency of the fast-axis direction of lithospheric anisotropy with large structures observed at the surface suggest a continuation of the main shear zones into the lithospheric mantle, as suggested by Vauchez et al. (2012). Nonetheless, a few exceptions are observed, for example, at stations sabr, sbbr and nbpb. Such discrepancies in the anisotropy orientations could be related to more local features such as fluid content, presence of cracks or plutonic bodies along the shear zones, fractures or mineral assemblages (Levin and Park, 1997; Mainprice and Nicolas, 1989).

5.2 Anisotropy along the passive margin

Inspection of stations located along the eastern and equatorial margins reveals that anisotropy exhibits – on average – directions either perpendicular or oblique to the coast in the lithospheric mantle. This observation is in agreement with SKS-splitting measurements in this area performed by Bastow et al. (2011) and Assumpção et al. (2011). These authors concluded that the anisotropy reported from SKS splitting along the northeastern Brazilian margins must be related to fossil anisotropy inherited from the opening of the South Atlantic Ocean. This interpretation is based on the relatively small time delay measured along the coast.

We compare the independent SKS-splitting measurements with our results from harmonic stripping of receiver functions in Fig. 7. For a better comparison we chose to represent the k=2 harmonics at stations where SKS splitting was measured because: (i) we expect only horizontal (recorded on the k=2 harmonics) or slightly dipping anisotropy in such geodynamical context; (ii) we observe in our data that k=1 and k=2 harmonics display energy within the same order of magnitude, suggesting slightly dipping to horizontal anisotropy beneath most stations; and (iii) SKS waves are mainly sensitive to (sub)horizontal anisotropy (Levin et al., 2007). Figure 7 shows a good agreement between anisotropic orientations recorded by receiver functions and SKS splitting along the eastern and equatorial margins, confirming that the recorded anisotropy beneath coastal stations is mainly located in the lithospheric mantle. The oblique to parallel orientation of anisotropy along the east and equatorial coasts, respectively, is consistent with the opening trend of the margin (Moulin et al., 2010).

5.3 Non-azimuthal anisotropy along the aborted Cariri–Potiguar rift

At a number of stations (ar05, nbma, pfbr, nbpa and cs6b), uncertainties regarding the direction of the fast axis of anisotropy are larger than 20^∘. These stations, however, display similar energy than stations with smaller uncertainties (see Fig. 6 and Fig. S3 in the Supplement). Interestingly, those stations seem to form a remarkable line trending northeast–southwest that approximately coincides with the location of the Cariri–Potiguar trend. Stations nbta and pcse also seem to align along the same direction more to the east.

This northeast–southwest oriented line is located above a northeast–southwest trending channel of thin lithosphere imaged by the tomographic study of Simões Neto et al. (2019). We suggest that deformation from thermomechanical erosion by horizontal, sublithospheric flow along the channel – also postulated by Simões Neto et al. (2019) – must be ongoing above this northeast–southwest channel. Furthermore, as initial thinning of the lithosphere along the channel was triggered by Mesozoic extension along the Cariri–Potiguar trend, alterations to the original Precambrian anisotropic fabric by Mesozoic extension might still be present. Additionally, we note that the location of the Cariri–Potiguar trend also marks the boundary between the east–west striking shear zones in the southern province from the northeast–southwest striking shear zones in the western province (Fig. 1). This suggests the Cariri–Potiguar trend also marks the location of a former paleo-suture that later acted as a zone of weakness along which the Mesozoic rift (now aborted) could develop. Thus, we believe the non-azimuthal anisotropy recorded at stations located along this trend is likely related to complex fossil anisotropic fabrics resulting from a combination of deformation along the ancient collision between Precambrian blocks, Mesozoic extension and thermomechanical erosion/mantle dragging by sublithospheric flow.

6 Conclusions

We have investigated depth-dependent anisotropy in the Borborema Province of northeast Brazil through harmonic decomposition of receiver functions developed at 39 stations in the region. Our main results include the following: (i) anisotropy within the province is characterized by a horizontal to slightly dipping fast axis of symmetry; (ii) consistency of anisotropic orientations within the crust and the lithospheric mantle suggest a continuation of surface shear-zones down to lithospheric depths; (iii) fast axes of symmetry are oriented parallel to the main shear zones within the continental interior and subparallel to Mesozoic extension along the passive margins, consistent with a fossil origin inherited from the opening of the South Atlantic Ocean; (iv) large uncertainties in the anisotropy orientation along a northeast–southwest trending line in the center of the province might be related to complex fossil anisotropic fabrics resulting from a combination of deformation along the ancient collision between Precambrian blocks, Mesozoic extension and thermomechanical erosion/mantle dragging by sublithospheric flow identified in an independent tomography study.