Lithospheric and sub-lithospheric deformation under the Borborema Province of NE Brazil from receiver function harmonic stripping

The depth-dependent anisotropic structure of the lithosphere under the Borborema Province of northeast Brazil has been investigated through 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 tec5 tonic processes. Our results reveal the presence of anisotropy within the crust and the lithospheric mantle throughout the entire Province, with the exception of a few stations in the continental interior that lack evidence for any anisotropic signatures. Most stations in the continental interior report consistent anisotropic orientations in the crust and lithospheric mantle, suggesting a dominant NE-SW pervasive deformation along lithospheric-scale shear zones developed during the Brasiliano-Pan African orogeny. The lack of anisotropy at a few stations along a NE-SW trend in the center on the Province is harder to explain, 10 but might be related to heating of the lithosphere by an asthenospheric channel. 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. 15 Copyright statement.


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 of NE Brazil, for instance, has witnessed several cycles of deformation, as well as recurrent episodes of intraplate volcanism and uplift, during its geological 20 history. Brasiliano-Pan African deformation is well represented through the network of shear zones that pervade the Borborema Province Neves et al., 2000). These shear zones separate several tectonic terrains of Paleoproterozoic and Archean age 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). 5 Also, current topography of the Borborema Plateau and the Sertaneja Depression may have resulted from a combination of on-going 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, b), ambient noise (Dias et al., 2014) or P-wave tomography (Simões Neto et al., 2019), and SKS splitting (Bastow et al., 2015) 10 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, figure 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 15 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 Mezosoic 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 20 mantle. Our results also confirm the absence of lithospheric (fossil) anisotropy under stations that did not record SKS-splitting.
These stations are aligned along a NE-SW trend located west and north of the Borborema Plateau, a high-standing topographic feature that rises ∼1000 m above mean sea level. We argue that the absence of anisotropy in the lithosphere is related to sub-lihospheric heating of the overlying lithosphere by a shallow asthenospheric channel under the Province.
2 Geological setting 25 The Borborema Province formed during the Neoproterozoic Braziliano/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). It thus 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 30 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 EW and NE-SW ( Figure 1). These shear zones are major structures several hundreds of kilometers long and tens of kilometers wide (Vauchez and da Silva, 1992) 2 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 . The shear zone network can be split into two domains: a western domain of NE-striking skike-slip faults, and an eastern domain of more sinuous, discontinuous E-W-striking shear zones . 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 NE 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 been traditionally proposed. On one hand, the accretionary model proposes that the Borborema Province comprises of several Paleoproterozoic small continental fragments that aggregated along the shear zones, which then constitute lithospheric-5 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 (Figure 1). 10 On the other hand, it has been suggested that the Borborema Province was 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.
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 30 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. On one hand, 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 35 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., 2015b, a). However, Luz et al. (2015b) debated the time of emplacement of such a 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 5 lineament postulated from seismic P-wave tomography (Simões Neto et al., 2019). The tomographic study also identifies an asthenospheric low-velocity channel trending NE-SW 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 NE Brazil. Hz. Further detail on the 75 seismic stations is given in the Table S1 and their geographical location displayed in Figure 2.

Receiver function processing and migration
Receiver functions were computed for the 75 stations making 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 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 5 6 functions 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-wavefront with crustal and upper mantle discontinuities under the recording station (Langston, 1979(Langston, , 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 5 (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.
The main processing steps involved in the development of the receiver function estimates are summarized below, and further details about 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 10 Figure 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 re-sampled 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 through the iterative, time-domain 15 procedure of Ligorria and Ammon (1999), with 500 iterations. The deconvolved time series were again filtered with the same Gaussian filter of 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 20 of the migration is to correct the phase move-out introduced by varying incidence angles among the incoming teleseismic Pwavefronts, effectively equalizing the receiver function waveforms in the depth domain (Dueker and Sheehan, 1997). Migration before harmonic stripping at individual stations was previously utilized by Audet (2015) decomposition on depth-migrated cross-sections obtained through CCP stacking of receiver functions. Next, the migrated radial 25 and transverse receiver functions for each station were grouped by back-azimuth in 36 non-overlapping, 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 9 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 9 bins (each 10°wide) allows the mapping of either half the period for a 2-lobed pattern (anisotropy with plunging fast axis of symmetry) or a full period for a 4-lobed 30 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 Figure 4.

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);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 5 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 +π/2k and naturally improve the azimuthal coverage around the station. After the harmonic decomposition is performed, up to 5 coefficient functions, corresponding to the first three harmonics, are obtained (k = 0, 1, 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 10 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 N-S and E-W 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 N-S 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 5 coefficients for the 3 harmonic degrees (k = 0, 1, 2) through a singular value decomposition. These harmonics are calculated for every depth within a selected depth-window. 5 The matrix equation that implements the harmonic decomposition is given by 0 cos(φ 1 + π/2) sin(φ 1 + π/2) cos(2(φ 1 + π/4)) sin(2(φ 1 + π/4)) . . . . . . . . . . . . . . . 0 cos(φ n + π/2) sin(φ n + π/2) cos(2(φ n + π/4)) sin(2(φ n + π/4)) 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 equation (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 5 displays non-zero 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): and These energy functions allow to discriminate between dipping interfaces/a plunging axis of symmetry and horizontal anisotropy. If E k=1 > E k=2 , the dominant anisotropy is either a dipping interface or an anisotropic layer with a plunging 10 axis of symmetry. In that case, we rotate B(z) and C(z) in discrete back-azimuth increments α (where α ∈ [0, 2π]) 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=2 > E k=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 α ∈ [0, 2π]) and search for the value of α that maximizes D(z) 15 (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.
An example of harmonic decomposition is shown for station PFBR in Figure 5. In order to estimate uncertainties, we applied a bootstrap statistical approach by randomly re-sampling with replacement our receiver functions. We performed such analysis with 200 replications at each of the selected stations. From these 200 values, we estimated the standard error (2σ), which 20 corresponds to the uncertainty in the direction of the fast-axis of symmetry. A measurement is considered as not reliable, and then rejected, if the estimated uncertainties are larger than 20°.

Results
Anisotropy parameters were examined for each station at two depth-window ranges: (1) crust ( Figure 6A), which was assumed to be located between 0 and 33 km depth, in agreement with the 32-40 km range estimated by Luz et al. (2015b) under the 25 Borborema Plateau and 30-33 km under the surrounding basins; and (2) lithospheric mantle, which was taken to be between 33 and 100 km depth ( Figure 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.
An inspection of Figure 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 section 5. Unresolved anisotropic directions within the crust are recorded around 5 longitude -40°for stations nbpb, ar02, ar05, ar06, nbpn, at the border of the Borborema Plateau (stations nbta, pctv, nbli, caub), and within the Sergipe-Alagoas and Pernambuco basins (stations nban and pcal). The majority of stations that sample clear anisotropic directions display a NE-SW to E-W trending axis of symmetry, except stations cs6b (trend ∼ NNW-SSE), km60 and nbma (trends N-S). 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, 10 pcsl, pcja, lp06, 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 supplementary materials, Figure S1). For example, station PFBR ( Figure 5) shows clear, non-zero 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 section 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), 5 and along a NE-SW axis located northwest of the Borborema Plateau (nbma, pfbr, nbpa, cs6b). Most anisotropic directions trend NE-SW to E-W, with the exception of stations km60 and nbma that show N-S 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, rcbr). Note that stations located along the continental margin show anisotropy within the lithospheric mantle with a fast axis of symmetry 10 that is oblique (stations nbmo, nbpv, nbrf, nban, nbit) or perpendicular (stations nbcl, 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).

Pervasive anisotropy with (sub)horizontal fast axis of symmetry
As described is section 4, we observe in northeast Brazil a dominance of 2π-periodicity (k=1) anisotropy in the lithospheric mantle, 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. 10 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 supplementary materials, Figure 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 sub-parallel to the main E-W to NE-SW shear zone directions 5 (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). A few exceptions, nonetheless, 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). 10

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 15 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 Figure 7. For a better comparison we chose to represent the k=2 harmonics at stations where SKS splitting were measured because: (i) we expect only horizontal (recorded on the k=2 harmonics) or slightly dipping anisotropy in such geodynamical 20 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 equa-25 torial coasts, respectively, is consistent with the opening trend of the margin (Moulin et al., 2010).

Non-azimuthal anisotropy along the aborted Cariri-Potiguar rift
At a number of stations (ar05, nbma, pfbr, nbpa, cs6b), uncertainties for the direction of the fast axis of anisotropy are larger than 20°. These stations, however, display similar energy than stations with smaller uncertainties (see Figure 6 and supplementary material Figure S3). Interestingly, those stations seem to form a remarkable line trending NE-SW that approximately coincides 30 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 NE-SW oriented line is located above a NE-SW trending channel of thin lithosphere imaged by the tomographic study of Simões Neto et al. (2019). We suggest that deformation from thermo-mechanical erosion by horizontal, sub-lithospheric flow along the channel -also postulated by Simões Neto et al. (2019) -must be ongoing above this NE-SW channel. Also, 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 EW striking shear zones in the southern Province from the NE-SW striking shear zones in the western Province ( Figure 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 5 related to complex fossil anisotropic fabrics resulting from a combination of deformation along the ancient collision between Precambrian blocks, Mesozoic extension, and thermo-mechanical erosion/mantle dragging by sub-lithospheric flow.

Conclusions
We have investigated depth-dependent anisotropy in the Borborema Province of NE Brazil through harmonic decomposition of receiver functions developed at 39 stations in the region. Our main results include: (i) anisotropy within the Province is