Crustal structure of the East African Limpopo margin, a strike-slip rifted corridor along the continental Mozambique Coastal Plain and North Natal Valley

Coincident wide-angle and multi-channel seismic data acquired within the scope of the PAMELA Moz3-5 project allow us to reconsider the formation mechanism of East African margins offshore of southern Mozambique. This study specifically focuses on the sedimentary and deep-crustal architecture of the Limpopo margin (LM) that fringes the eastern edge of the Mozambique’s Coastal Plain (MCP) and its offshore southern prolongation the North Natal Valley (NNV). It relies primarily on the MZ3 profile that runs obliquely from the northeastern NNV towards the Mozambique basin (MB) with additional inputs from a tectonostratigraphy analysis of industrial onshore–offshore seismic lines and nearby or crossing velocity models from companion studies. Over its entire N–S extension the LM appears segmented into (1) a western domain that shows the progressive eastward crustal thinning and termination of the MCP/NNV continental crust and its overlying pre-Neocomian volcano-sedimentary basement and (2) a central corridor of anomalous crust bounded to the east by the Mozambique fracture zone (MFZ) and the oceanic crust of the MB. A prominent basement high marks the boundary between these two domains. Its development was most probably controlled by a steep and deeply rooted fault, i.e., the Limpopo fault. We infer that strike-slip or slightly transtensional rifting occurred along the LM and was accommodated along this Limpopo fault. At depth we propose that ductile shearing was responsible for the thinning of the continental crust and an oceanward flow of lower crustal material. This process was accompanied by intense magmatism that extruded to form the volcanic basement and gave the corridor its peculiar structure and mixed nature. The whole region remained at a relative high level during the rifting period and a shallow marine environment dominated the pre-Neocomian period during the early phase of continent–ocean interaction. It is only some time after break-up in the MB and the initiation of the MFZ that decoupling occurred between the MCP/NNV and the corridor, allowing for the latter to subside and become covered by deep marine sediments. A scenario for the early evolution and formation of the LM is proposed taking into account both recent kinematic and geological constraints. It implies that no or little change in extensional direction occurred between the intra-continental rifting and subsequent phase of continent–ocean interaction.

1 Introduction

The break up of East Gondwana led to the separation of four independent continental blocks, namely Africa, Antarctica–Australia, Madagascar–India and Patagonia. This fragmentation strongly segmented the East African margins into a succession of divergent and strike-slip segments (Fig. 1). To the north, southward motion of the Madagascar–India plate along the Davie Fracture Zone (DFZ) opened the Somali Basin (SB). Similarly, the Antarctica–Australia plate drifted southward with respect to Africa along the Mozambique Fracture Zone (MFZ) and opened the Mozambique Basin (MB). Finally, further southwest, the Patagonia plate escaped from the coast of South Africa along the Agulhas-Falklands Fracture Zone (AFFZ) while seafloor spreading occurred in the South Natal Valley (SNV).

Figure 1 Free-air satellite-derived gravity anomaly of the Mozambique Channel (Sandwell et al., 2014) showing the complexity and segmentation of the East African margin (modified from Thompson et al., 2019). Segments 1 and 2 correspond to the movement of the Madagascar and Antarctica tectonic plates that opened the Somali Basin (SB) and the Mozambique Basin (MB), respectively. Segment 3 corresponds to the drifting of the Falkland and Patagonian plates that opened the South Natal Valley (SNV). NNV stands for North Natal Valley, MCP stands for Mozambique Coastal Plain, NMR stands for North Mozambique Ridge, SMR stands for South Mozambique Ridge, NMadR stands for North Madagascar Ridge and SMadR stands for South Madagascar Ridge. Our study area is the Limpopo margin (LM, in bold). In red is volcanism without distinction about the age, black lines are major tectonic features and brown lines denote craton contours (Mercator projection).

The exact timing and mechanisms of Gondwana breakup remain, however, still speculative. Namely, the Mozambique Coastal Plain (MCP) in southern Mozambique and its offshore prolongation the North Natal Valley (NNV) form a wide buffer zone inherited from this complex breakup (Figs. 1 and 2). It is surrounded to the north and west by the Lebombo and Mateke–Sabi monoclines, representing continental flood basalts (CFBs) and associated dike swarms of Lower Jurassic Karoo age (e.g., Jourdan et al., 2005; Watkeys, 2002), and to the south and east by oceanic crust of the MB and SNV of Upper Jurassic (e.g., Mueller and Jokat, 2019) and Lower Cretaceous (Goodlad et al., 1982) ages, respectively. In the absence of modern deep seismic investigation, the crustal nature of the area is highly debated. Many have interpreted the presence of thickened oceanic crust and/or thinned continental crust with high magmatic content based on scarce seismic data, potential field data, and/or geological correlation with the conjugate margin in Antarctica’s western Dronning Maud Land (DML) (e.g., Klausen, 2009; Leinweber and Jokat, 2011; Tikku et al., 2002; Watts, 2001). Consequently, plate kinematic frameworks including so-called “tight-fit” reconstruction have been privileged, where the DML fully or partly overlaps the MCP and NNV and an initial phase of rifting oblique or normal to the subsequent southward plate drift is inferred (e.g., Martin and Hartnady, 1986; Cox, 1992; Klausen, 2009; Mueller and Jokat, 2019).

To shed new light on the East African margin structures and their formation processes, the PAMELA (PAssive Margins Exploration Laboratories) project was initiated in the early 2010s by TOTAL and IFREMER in collaboration with French universities (Université de Bretagne Occidentale, Université Rennes 1, Université Pierre and Marie Curie), the CNRS and the IFPEN. In 2016, the PAMELA MOZ3-5 cruise acquired seven combined multi-channels and wide-angle seismic profiles in the MCP and NNV area (Moulin and Aslanian, 2016; Moulin and Evain, 2016) (Fig. 2). Velocity models built from these data challenge the tight-fit scenario as they unravel the continuity of a thick and highly intruded continental crust beneath the MCP and NNV (Moulin et al., 2020; Leprêtre et al., 2021; Watremez et al., 2021). Conversely, these results strengthen previous studies that interpreted continental crust in the area (Domingues et al., 2016; Hanyu et al., 2017) and brings strong support for a “loose” fit plate reconstruction without plate overlap (Thompson et al., 2019). Within such a framework, both rifting and subsequent spreading are inferred to proceed from a continuous quasi N–S movement of the Antarctica–Australia and Madagascar–India plates with respect to Africa with major implications for the dynamics of rifting along the East African margins and the age of Gondwana breakup.

Figure 2 Maps of the PAMELA MOZ35 seismic experiment. Combined MCS and wide-angle profiles are shown as black lines. MZ3 is highlighted, with red dots marking positions of OBS. Instruments (OBS and land stations) and names of other profiles are in white. (a) Topography and bathymetric data (Smith and Sandwell, 1997) with inset for location of studied area. The gray stars indicate the location points of 1D Vs depth profiles from Domingues et al. (2016). LM is the Limpopo margin, shown in red; NNV is the North Natal Valley; N-SNV is the northern South Natal Valley; and MR is the Mozambique Ridge. (b) Free-air gravity anomaly grid (Pavlis et al., 2012; Sandwell et al., 2014). Main features known in the area are from Mueller and Jokat (2019), which are based on Leinweber and Jokat (2012) and Mueller and Jokat (2017) for the magnetic spreading anomalies in the Mozambique Basin, and on Goodlad et al. (1982) for the magnetic anomalies in the South Natal Valley. The outline of the Naude Ridge (NR) is from Goodlad (1986). ALR is the Almirante Leite Ridge, AG is the Ariel Graben, DP is the Dana Plateau, GP is the Galathea Plateau, and VP is the Vauban Plateau.

To further constrain the Gondwana breakup, this study focuses on the MZ3 profile acquired during the PAMELA MOZ3-5 experiment (Fig. 2). It investigates how rifting occurred to the east of the MCP/NNV along the Limpopo margin (LM), one of the East African strike-slip segments. We produce both a P-wave velocity model of the sedimentary and deep crustal structure along the profile and a tectonostratigraphic analysis of the area based on key industrial seismic lines. Our results are interpreted jointly with models from companion profiles to reveal the peculiarities of this margin and propose a scenario for its formation and early evolution within tight kinematic constraints.

2 Geological setting and controversies

East Gondwana breakup initiated in a heterogeneous lithosphere that was last reworked and assembled during the Pan-African Orogeny at about 720 to 550 Ma (Guiraud et al., 2005; Jacobs and Thomas, 2004). In southern Africa, Pan-African structures include the Zambezi, Lurio and Mozambique belts. They linked together Archean and Proterozoic cratons (Zimbabwe and Kaapvaal) and belts (Limpopo, Namaqua-Natal) to their equivalents in Antarctica's western Dronning Maud Land, namely the Grunehogna Craton and Maud Belt (Bingen et al., 2009; Jacobs et al., 2008; Riedel et al., 2013).

The Karoo sedimentary sequence and subsequent Karoo magmatism attest to widespread intercontinental tectonic and magmatic activities before Jurassic break-up. Karoo sediments consists of Late Paleozoic to Early Jurassic deposits within marginal basins in South Africa (Catuneanu, 2004) and interconnected rifts in Antarctica, eastern Madagascar, and central and East Africa (Daly et al., 1989; Elliot and Fleming, 2004; Geiger et al., 2004; Salman and Abdula, 1995). They were largely covered by Karoo continental flood basalts between 185 and 177 Ma while up to 172 Ma sills and dikes also massively intruded Precambrian basement structures in southern Africa and Antarctica's Dronning Maud Land (Cox, 1992; Hastie et al., 2014; Jourdan et al., 2007).

The MCP is known from industrial wells to be floored by volcanics but solid dating is lacking in support for a suspected Karoo age (Flores, 1973; Salman and Abdula, 1995). The oldest reliable age for basalts in the area concerns the Lower Cretaceous Movene igneous that disconformably cover Karoo sediments to the east of the Lebombos monocline (Flores, 1973; Watkeys, 2002). As already stressed in companion papers (e.g., Moulin et al., 2020; Leprêtre et al., 2021), the crustal nature in the area is the subject of a vigorous debate. The MCP gravity signature was previously interpreted as suggestive of oceanic crust (Watts, 2001), while continental crust was more recently favored based on ambient noise tomography (Domingues et al., 2016). Controversies also exist regarding the interpretation of weak magnetic anomalies within the NNV, which could either sign fossil spreading centers (Green, 1972; Tikku et al., 2002; Martin and Hartnady, 1986; Leinweber and Jokat, 2011; Mueller and Jokat, 2019) or magmatic intrusions within a thinned continental crust (Hanyu et al., 2017). Available marine seismic data in the NNV and its conjugate margin were used in the past in support of either continental crust (Dingle and Scrutton, 1974; Lafourcade, 1984) or oceanic crust (Ludwig et al., 1968; Jokat et al., 2004; Green, 1972). With new deep seismic data acquired within the scope of the PAMELA project, companion papers of this study argue for continental crust below the MCP and NNV. Despite its high magmatic content, Moulin et al. (2020) and Leprêtre et al. (2021) show the continuity of a thick crust with a comparable velocity structure to surrounding cratons and other continental structures worldwide.

Controversies further exist regarding the age of the first oceanic crust within the MB along the eastern MCP and NNV. It is generally considered to be around 155 Ma on both the Mozambique margin and its conjugate Riiser-Larsen Sea in Antarctica (Jokat et al., 2003; König and Jokat, 2010; Thompson et al., 2019). However, some studies reviewed the oldest magnetic anomalies found on the northeastern corner of the MB together with evidence of fracture zone paths on the gravity field to postulate an earlier initiation of seafloor spreading at 166–164 Ma (Leinweber and Jokat, 2012; Mueller and Jokat, 2019). This situation slightly differs on the opposite corner of the MB, where early rifting concentrated in the offshore Zambezi depression before a rift jump isolated the Beira continental block from the rest of margin (Mahanjane, 2012; Mueller et al., 2016). Just south of the Beira High (BH), the oldest magnetic anomaly identified is therefore dated at 155–157 Ma (Mueller and Jokat, 2019) (Fig. 2).

In between the MCP and NNV complex and the MB, the interpretation of a clear N–S trending positive gravity anomaly is disputed. For many it is the natural prolongation of the MR oceanic plateaus located to the south (Fig. 2). However, the tectonostratigraphic study of Li et al. (2021) suggests a volcano-sedimentary basement here that formed during the intra-continental rifting period. In this case, the MR would not extend further north, but instead positive gravity anomalies would reflect the presence of large sedimentary features such as contourites, which are commonly observed in the area (Moulin et al., 2020; Thiéblemont et al., 2019).

All of these elements have major implications regarding the location of the continent–ocean transition (COT). For Moulin et al. (2020), the SW–NE-oriented Ariel Graben (AG) and Naude Ridge (NR) mark this COT. In tight-fit kinematic models, however, it is located just east and south of the Lebombo and Mateke–Sabi Monoclines, respectively, or slightly further inside the MCP (Martin and Hartnady, 1986; Klausen, 2009; Leinweber and Jokat, 2011; Mueller and Jokat, 2019). According to Mueller and Jokat (2019), the COT could also be located in the northern part of the NNV where several southwards ridge jumps may have occurred. Alternatively, Hanyu et al. (2017) divide the MR into a northern continental part, including the Dana and Galathea plateaus, and a southern oceanic part, including the Vauban plateau (Fig. 2). This places the COT even further south than the AG and NR despite several studies that suggest the crustal structure and morphological characteristics of these plateaus formed in Lower Cretaceous on a triple junction (Fischer et al., 2017; Gohl et al., 2011; König and Jokat, 2010; Leinweber and Jokat, 2012). It should be noticed, however, that some dredges made along the steep eastern border of the MR facing the MFZ and on its southwestern edge facing the SNV recovered Archean fragments of continental rocks and metamorphic samples of African affinity (Ben-Avraham, et al., 1995; Hartnady et al., 1992; Mougenot et al., 1986).

Following the formation of the southern Mozambique margins and MR, intense post-rift tectonic activities reportedly affected the area. Long-term subsidence is attested by the thick cover of Cretaceous to Neogene sediments over the MCP and NNV, although this was perturbed by several tectono-volcanic events (Baby, 2017). Cenozoic to modern extension has also been emphasized within the MCP (Mougenot et al., 1986) and possibly southward in the NNV (Wiles et al., 2014). These days extensional stress fields are reported in both southern Africa and Madagascar (Adam and Lebedev, 2012; Piqué et al., 1999; Tsang-Hin-Sun et al., 2021), while there is increasing evidence for the southward propagation of the East African rift system through both onshore and offshore branches (Mougenot et al., 1986; Deville et al., 2018; Courgeon et al., 2018). Those branches where recent volcanism has concentrated possibly delineate actual micro-plate boundaries (e.g., Saria et al., 2014).

3 Combined multichannel seismic (MCS) and wide-angle seismic (WAS) acquisition, processing and modeling

3.1 Generalities

In 2016, the PAMELA MOZ3-5 cruises onboard French R/V Pourquoi Pas? acquired a total of seven coincident multi-channel and wide-angle seismic profiles over the southeastern MCP, the offshore NNV and the LM (Fig. 2; Moulin et al., 2020; Moulin and Aslanian, 2016; Moulin and Evain, 2016). This study concentrates on the 500 km long MZ3 profile trending SW–NE offshore of the MCP. It runs from the NE edge of the NNV and crosses the LM towards the northwestern corner of the MB. To illuminate the deep crustal structure, 32 four-component ocean bottom seismometers (OBSs) from Ifremer’s Marine Geosciences Unit were deployed every 7 NM. A total of 2960 shots were recorded by those instruments and by a 4.5 km, 720-channel towed marine streamer. MZ3 crosses four out of the seven profiles acquired during the MOZ3-5 cruise: MZ2 and MZ6 (Schnürle et al., 2018) acquired over the NNV and MZ4 and MZ5 that also cross the LM (Watremez et al., 2021). Two other profiles acquired in the NNV, MZ1 and MZ7, are presented by Moulin et al. (2020) and Leprêtre et al. (2021).

Ifremer’s SolidQC software was initially used for processing MCS data. It allows for data quality control, 2D geometry setup and SEG-Y file generation. Further processing was performed using CGG-Veritas Geocluster software. This sequence included external mute of direct and water wave arrivals, large band-pass filtering (1, 8, 64 and 92 Hz edge frequencies), predictive deconvolution (440 ms operator, 252 ms IO, 32 ms gap), 4 ms resampling, seafloor multiple attenuation using surface-related multiple elimination (SRME) and Radon transform (RAMUR), time-variant band-pass filtering, and pre-stack Kirchhoff migration. Velocity picking is done after each major processing step in order to refine the final velocity model and build a coherent pre-stack Kirchhoff time-migrated section (PSTM).

OBS WAS data processing included clock drift corrections, instrument localization using direct water wave arrivals to correct for the drift from their deployment position and eventually band-pass Butterworth filtering to improve travel time arrival identification and picking.

We used the iterative two-dimensional forward ray-tracing procedure followed by the damped least-squares travel time inversion of the RAYINVR software (Zelt and Smith, 1992). Our modeling proceeded following a top-to-bottom strategy of arrival times fitting of both wide-angle reflected and refracted phases identified on OBS data. Starting from an initial velocity model including the bathymetry, model interfaces were inserted when a velocity change in apparent refracted velocity and a coincident high-amplitude reflective phase were observed (i.e., phase triplication). Layer velocities and thicknesses were adjusted in order to fit these phases. By converting our velocity model from depth to two-way travel time we checked and adjusted the geometry of the upper sedimentary and acoustic basement interfaces as identified on the MCS section (see Fig. 3b).

Below we describe in detail what guided our modeling by describing layers and their facies observed on the MCS line, followed by their corresponding refracted and reflected phases identified in OBS records. Overall, by combining MCS and WAS data we identified four geological units along MZ3. From top to bottom they consist of (a) the upper sedimentary package (orange and yellow Ps phases or interfaces in Figs. 3 and 9); (b) the volcano-sedimentary sequence that forms the acoustic basement on the MCS (green Psv phases or interfaces); (c) the crystalline basement (blue Pg phases or interfaces); and (d) the upper mantle (magenta Pn phases or interfaces). From MZ3 WAS data, this represents a total of 72 750 picked arrival times that are summarized in Appendix Tables A1. Travel time uncertainty for each phase was automatically calculated based on the trace signal-to-noise ratio. They range from 0.025 s for high ratios to 0.25 s for poor ratios.

Figure 3 (a) Two-way travel time record section of MCS data along the MZ3 profile. (b) The same as (a) but overlain by time-converted interfaces from the coincident wide-angle velocity model. Intersections with other profiles from the MOZ35 cruise are indicated by a red line. OBS locations are indicated by inverted triangles.

3.2 MZ3 MCS line

MZ3 PSTM is shown in Fig. 3. The first 100 kilometers of the MZ3 profile run over the NNV are dominated by volcanic edifices of the Almirante Leite Ridge (ALR) which completely masks the underlying structure (−20 and 40 km). On each side of the blind zone, a maximum of 1 s twt (second two-way travel time) of upper sediments is present. These consist of a first layer (S1) of weak signal and poorly resolved horizons below the seafloor. It is very thin to the SW of the ALR but thickens on its opposite side. At the bottom of S1 a stronger reflector (top-S3) is flat and parallel to the seafloor to the SW of the ALR but dips and undulates towards the NE. A few strong amplitudes and unconformable reflections exist on each side of the ALR beneath this interface. A set of more continuous and coherent horizons marks the top of acoustic basement (top-SV1), which appears highly and steeply fractured to the SW of the ALR. At least 0.5 s twt of parallel layering can be identified within the upper basement (SV1) on this side, while a more chaotic and incoherent signal dominates on the NE side of the ALR. Because energy is lost below about 3 s twt, the SV2 and deeper layers are constrained by OBS data.

With increasing model distance, the seafloor over the LM forms a slight dome between 120 and 260 km and then continuously deepens. S1 continues from the NNV as a 0.5 s twt thick layer, the internal structure of which is almost transparent because it is poorly resolved by our seismic acquisition. Below, S2 is a layer of intense but highly disrupted reflectivity revealing various sedimentary features, possibly including channel incisions (100 km model distance), contourites (360–400 km) and mass transport deposits (400–430 km). This contrasts with the S3 layer, where more continuous and organized horizons are present, though these are of weaker amplitude. The top of the acoustic basement (top SV1) corresponds to a set of parallel and stronger-amplitude reflectors in continuity with the NNV up to 300 km model distance. There, a transition to a rugged surface marked by a broad basement high occurs. To the SW of the high, the basement internal structure (SV1) shows clear parallel or dipping reflectors, while to its NE it is chaotic with incoherent signal. There are two trends in dipping reflectors: one to the NE that is slowly dipping toward the western edge of the basement high (240–260 km) and a second, dipping to the SW, that coincides with highly fractured and sharp basement. Only the internal structure of the SV1 layer, about 0.5 s twt thick, is imaged by the MCS. No energy reflects further down, and SV2, the crystalline basement (G1–G4), and the Moho (top-Moho in Fig. 3b) are constrained by wide-angle seismic acquisition only.

3.3 MZ3 WAS data

We selected key OBS sections to illustrate this section, but all records can be found in Supplement with their corresponding synthetic sections, travel time picks and fits, and ray-tracing through the final model and the coincident portion of MCS line.

Figure 4 Sedimentary and volcano-sedimentary seismic phases of MZ3OBS03 (top) and MZ3OBS14 (bottom). For each instrument, the six panels display the following data: (a) seismic record, (b) seismic record with color-coded predicted arrivals, (c) synthetic section with color-coded predicted arrivals, (d) color-coded picked travel times with uncertainty bars overlain by dotted predicted times, (e) color-coded seismic rays, and (f) MCS time migrated section with color-coded model interfaces. (a–d) Travel time is reduced by a velocity of 5 km s⁻¹.

On most OBS records, two or three refracted and associated reflected phases describe the sedimentary package (e.g., orange phases in Figs. 4–5). Those refracted phases have apparent velocities and offsets rarely exceeding 3 km s⁻¹ and 20 km, respectively. They are visible as first arrivals at very short offsets then generally extend as secondary arrivals. Corresponding reflected phases are generally well observed at short offsets inside the cone formed by direct wave arrivals (red phase). Their travel times are coherent with those of the main horizons seen on the MCS profile (top-S2 and top-S3 in Fig. 3b). Clear examples of sedimentary phases can be observed on MZ3OBS14 and MZ3OBS30 record sections (Figs. 4b and 5b) located on the LM. Sedimentary phases recorded on instruments located on top or near the ALR are highly heterogeneous. The velocity structure is difficult to model and thus poorly constrained. As shown on MZ3OBS03 (Fig. 4a) there are high apparent velocities at short offsets followed by a shadow area. This is modeled by a negative velocity contrast reflecting the presence of high-velocity extruded volcanic material on top of lower-velocity sediments.

Figure 5 The same as Fig. 4. Sedimentary and volcano-sedimentary seismic phases in MZ3OBS20 (top) and MZ3OBS30 (bottom).

The top of the acoustic basement is marked on OBS records by a strong-amplitude reflected phase (green Top-SV1), generally the latest visible within the water cone which correlates well with observations made on MCS (Fig. 3b). The amplitude contrast is created by a strong velocity jump of at least 1.5 km s⁻¹. Instruments located over the NNV and the LM west of the basement high indicates the presence of two layers within the acoustic basement. There is a first set of associated refracted and reflected phases (SV1) with apparent velocities lower than 5 km s⁻¹ and a second set (SV2) with velocities higher than 5 km s⁻¹. (e.g., MZ3OBS20 in Fig. 5a). The first of these layers when converted to two-way time coincides with the deep layering previously described in MCS (SV1 layer in Fig. 3b). However, its base and the SV2 layer below are entirely constrained by wide-angle data. To the NE of the basement high, OBS shows a clear change in the characteristics of the acoustic basement internal structure, as already suggested from its facies on MCS. Here a single layer (SV1) continues, while SV2 is absent below the basement high. Apparent velocities are between 5 and 5.5 km s⁻¹ (e.g., MZ3OBS30 in Fig. 5b), which are intermediate between those of SV1 and SV2 beneath the NNV and LM. Therefore, despite its prolongation toward the northeastern end of the line, the basement high appears as a clear marker of a change in nature of the acoustic basement.

Figure 6 Seismic phases from the crust and mantle on MZ3OBS25. The left-hand panels represent negative offsets (toward the SW), and the right-hand panels represent positive offsets (toward the NE). On each side, the six panels display (a) seismic record, (b) seismic record with color-coded predicted arrivals, (c) synthetic section with color-coded predicted arrivals, (d) color-coded picked travel times with uncertainty bars overlain by dotted predicted times, (e) color-coded seismic rays, and (f) MCS time-migrated section with color-coded model interfaces. (a–d) Travel time is reduced by a velocity of 7 km s⁻¹.

The crystalline basement of the LM in the central portion of MZ3 is very well constrained by OBS records. They show long offsets refracted phases with only mild variations in apparent velocities (blue Pg phases). These are associated with numerous and quite high amplitude reflected phases (blue top-Pg phases). These combined observations suggest a thick crust characterized by a low-velocity gradient but a strong internal reflectivity. At least five layers can be identified based on an average account of triplications between refracted and reflected phases. The apparent velocity of refracted phases ranges from 6.5 km s⁻¹ for the shallow crust to 7 km s⁻¹ for slightly above for the deep crust (e.g., MZ3OBS11 and MZ3OBS25 in Figs. 6 and 7). Triplication between crustal and mantle phases occurs typically at about 150 km offset. This offset decreases to less than 100 km for instruments located close to the NE end of the profile (e.g., MZ3OBS31 in Fig. 8), implying important crustal thinning.

Figure 7 The same as Fig. 6 but for MZ3OBS11.

A different pattern of crustal phases is indeed observed on OBS records from the NE extremity of the LM. Crustal refracted phases at short offset suggest a stronger gradient and an average apparent velocity around 6.5 km s⁻¹. The gradient reduces at greater offsets but apparent velocities do not exceed 7 km s⁻¹. Some reflected phases can be distinguished, but overall there is no more evidence for intense crustal layering, the most prominent corresponding to the Moho. The transition to mantle phases occurs at a relatively short offset (about 50 km on MZ3OBS31; Fig. 8) suggesting thin crust. As for the acoustic basement, all these features combined depict a clear change in crustal nature.

Figure 8 The same as Fig. 6 but for MZ3OBS31.

The deep velocity structure beneath the NNV is only partly constrained because volcanism at the ALR creates a seismic blind zone. Moreover, our acquisition geometry does not allow for optimal ray coverage beneath the edifice, with OBS only located on top and to the NE of the ridge. Thus the two upper crustal layers of the model are only reasonably well resolved, with clear refracted and reflected arrivals of apparent velocities around 6 and 6.5 km s⁻¹. The deeper crustal structure is poorly constrained (see MZ3OBS11 in Fig. 7). Though refracted phases cannot be identified, some higher-amplitude phases can be distinguished at specific offset ranges. Those can be explained as resulting from reflected waves from a layered medium similar to the one modeled for the LM. We can further interpret some far offset phases (e.g., MZ3OBS10 in Fig. 9) as reflection from the Moho if considering a depth of about 35 km for this interface as suggested by the modeling of the crossing MZ2 and MZ6 profiles (Schnürle et al., 2018). Nonetheless, the poor constraints, we can reasonably consider the deep crustal structure beneath the NNV and ALR as similar to the western part of the LM up to the basement high. Our dataset reflects a thick crust with strong internal reflectivity, a low-velocity gradient and high bottom velocities for both domains (7.3 km s⁻¹).

Figure 9 The same as Fig. 6 but for MZ3OBS10.

Arrivals from the mantle are observed on numerous OBS records throughout the profile. With many reflected phases converging towards the triplication point between crustal and mantle phases due to low-velocity gradient within the crust, there are generally high-amplitude signals at this place on OBS records, making the identification of the starting point of the mantle-refracted phase difficult. The consequence might be some larger uncertainties on the depth of the Moho and uppermost mantle velocities. However, these can be balanced using constraints from crossing profile MZ2 (Schnürle et al., 2018), MZ4 and MZ5 (Watremez et al., 2021), which agree on the depth of the Moho being a maximum of ± 2 km. Refracted mantle phases show increasing apparent velocities from 8.0 to 8.3 towards the NE (Figs. 6–9). Travel time fit for these phases is not perfect and is uneven from one record to another, which might also suggest important mantle heterogeneity along the profile. This heterogeneity is further evidenced by the presence of high-amplitude mantle-reflected phases on some records. Their dips are difficult to model with only horizontal layering. They may reveal local highly reflective body-diffracting seismic waves. Regardless, this suggests that the heterogeneous and highly reflective character of the crystalline basement might continue downward within the mantle of NNV and western LM. Again, we do not have evidence for such mantle reflectivity below the thin crust in the NE extremity of the LM, but we lack of offset on these OBS records.

4 Model assessment

Our final model explains 89 % (64 900) of the events picked with a global root-mean-square (rms) travel time of 65 ms, which given the uncertainty assigned for each phase results in a normalized chi-squared value of 0.65. Travel time rms and chi-squared values for each phase and each instrument are detailed in Appendix Tables A1 and A2, respectively. This section aims to present various elements that were used in order to validate our modeling and the final velocity model. They include a resolution analysis, our model's gravity response and its estimated depth and velocity uncertainties. Finally, a pre-stack depth-migrated (PSDM) MCS section is also presented.

Figure 10 Evaluation of the wide-angle model MZ3. (a) Model parameters, including interface depth nodes (squares), top- and bottom-layer velocity nodes (red circles); interfaces where reflections have been observed on OBS data are highlighted in blue. (b) Hit count for velocity (colored) and depth nodes (squares). (c) Spread function (SPF) values for velocity (colored) and depth nodes (squares). (d) Resolution of velocity (colored) and depth nodes (squares). Zones that were not imaged are left blank.

4.1 Model evaluation and resolution

Figure 10 presents four indicators of the quality of the MZ3 velocity model based on wide-angle data only. Interface depth node spacing and velocity node spacing are indeed key to model the lateral variations of the seismic velocity with sufficient resolution without introducing complexity not required by the data. Note that the four indicators are not calculated for the upper sedimentary package or the top of SV1 as those layers are constrained by coincident MCS data. Their topography is directly sampled on the MCS profile while velocities are modeled from wide-angle phases.

Logically, the density of velocity and depth nodes and the number of reflective segments are higher in the central part of the model and decrease towards its edges and with depth (Fig. 10a). Indeed, these regions of the model are less sampled by rays (Fig. 10b), i.e., usually from a single direction, and data quality is either poorer or degrades naturally with increasing offsets. Because the model parameterization has been adapted to these limitations, MZ3 shows overall limited smearing (± 3, Fig. 10c) and very good resolution (above 0.9, Fig. 12d). The quality of the modeling decreases at depth and towards model extremities. Resolution values remain, however, higher than 0.5, which are still considered acceptable. Greater smearing occurs essentially in the lower crust where only few refracted rays travel. Those have been difficult to identify on OBS records because many phases converge around to the triplication point with the PmP and Pn (see Sect. 3.3).

Figure 11 Gravity modeling along the MZ3 profile. (a) Density model showing blocks of equal density converted from the final velocity model according to empirical function (Ludwig et al., 1970). (b) Gravity anomaly along MZ3 profile. The free-air gravity anomaly derived from satellite altimetry (Pavlis et al., 2012) along the profile is in red, while yellow curves correspond to similar profiles extracted at 10, 20 and 30 km laterally (yellow lines). The blue curve is the gravity anomaly measured during the MOZ35 experiment. Green curves are the calculated gravity anomaly from the converted model in (a) according to several mantle densities (green labels).

4.2 Gravity modeling

We tested the gravity response of our final model against the measured and satellite-derived free-air gravity anomaly along the profile (Fig. 11). A 2D model consisting of homogeneous density blocks was constructed from the MZ3 final velocity model by converting P-wave velocities to densities according to empirical function (Ludwig et al., 1970). The resulting density ranges from 1800 kg m³ to 2400 kg m³ in the upper sedimentary units, 2500–2600 kg m³ in the volcano-sedimentary sequence and 2700 to 3000 kg m³ in the crystalline basement, while we also tested mantle densities between 3100 to 3300 kg m³. The model was extended down to 100 km where isostatic compensation may be reached and 500 km laterally on each side to avoid edge effects.

The gravity anomaly is relatively flat along and parallel to the MZ3 profile and does not exceed ± 30 mGal (Fig. 11b). There is sharp positive peak at the ALR and a broad negative one beneath the basement high of the LM that reflects the change in crustal structure. Between these two peaks the lithosphere appears in isostatic equilibrium. The calculated gravity anomaly based on the conversion of MZ3 velocity model reflects the overall shape of the observed anomaly all along the profile. However, an acceptable fit can only be obtained by adjusting converted densities. This is achieved by lowering the density of the mantle layers from 3300 to 3150 kg m³. Reduced mantle densities may signify the presence of a hotter than normal mantle. Alternatively, this may indicate that the conversions of crustal velocities into densities are underestimated under the NNV and western LM and overestimated at the NE extremity of the model. Our model does not generate a broad negative anomaly at the transition beneath the basement high of the LM, but the change in crustal properties generates a clear jump from negative to positive values consistent with the observed gravity anomaly.

Figure 12 Global uncertainty plots for MZ3. (a) Maximum and (b) minimum admissible velocity deviations from the preferred model, built from 126 models within the thresholds defined in the text. Shaded areas indicate ray coverage. Preferred model interfaces are indicated by black lines.

4.3 Uncertainties estimation

Velocity and depth uncertainties of MZ3 final velocity model were estimated using the VMONTECARLO code (Loureiro et al., 2016). VMONTECARLO explores the model space by generating random models. It evaluates the ability of each model to fit the observed dataset and translate it to estimates of uncertainties given some quality thresholds.

Due to computational costs, the explored model space was reduced by limiting the number of parameters and fixing some bounds. For MZ3 we chose to maintain fixed sedimentary and volcano-sedimentary layers, as well as basement depth nodes, because those benefit from refined MCS constraints. This results in a total number of variable parameters of 148 as shown in Fig. B1. We also allowed a maximum fluctuation of ± 0.50 km s⁻¹ on velocity nodes while the maximum depth node variations were set to bands of 1, 2 and 3 km for the upper crustal, the 3 mid-crustal and the Moho interfaces, respectively. We further limited the search by generating a maximum of 50 000 random models and imposed a scaling factor for velocity and depth bounds that, starting at 20 % of their maximum value, progressively increased to 100 % within the first half of the randomly generated models. Finally, the total number of observed travel times was reduced to be less than 50 000 picks to reduce ray-tracing computing time.

Scores are calculated for each model according to a function that takes into account travel time fit and the ratio of traced rays to the number of observed events (Loureiro et al., 2016). The MZ3 final velocity model gives a score of 0.88 with 43 928 rays traced over the 49 725 picks, an rms of 0.07 s and (after adjusting the pick's uncertainties) a χ² of 1. Among the 50 000 randomly generated models within the bounds given above, 46 110 were valid and scored up to 0.81. The Supplement gives a summary of this model space by plotting every 20 km model distance all 1D velocity–depth profiles color coded according to their normalized average score. By setting a confidence threshold at 95 % of the maximum normalized average score, local depth uncertainties for a particular velocity can be estimated (and vice versa). VMONTECARLO analysis shows that velocity uncertainties are globally stable around 0.1 km s⁻¹ both along the model and with depth. Similarly, depth uncertainties are on average quite stable along the profile. They increase with velocity from 0.5 km at 6.5 km s⁻¹ to 1.5–2 km at 8 km s⁻¹. Higher uncertainty values are observed around 300 km model distance because of the sharp transition in crustal structure and the potential lateral effects due to the orientation of MZ3. It is interesting to note that depth uncertainties for a velocity of 7 km s⁻¹ are usually high over a couple of kilometers to a maximum of 5 km in a few places. Up to 300 km model distance this can be explained by the very low-velocity gradient required in the middle to lower crust. In fact, VMONTECARLO analysis further confirms that most of the basement is composed of material with velocities higher than 6.5 km s⁻¹ even towards the NE end of the profile (distance greater than 320 km).

In order to visualize uncertainty estimates on profile a subset of possible alternative models were selected. These models respect the following criteria that we judge acceptable: a score over 75 % of the preferred model’s score, a χ² lower than 2, an rms lower than 0.1 s and at least 80 % of the rays traced by the preferred model. This subset represents 126 models that were combined to produce the minimum and maximum admissible velocity deviations maps shown in Fig. 12. Large uncertainties within the depth bounds allowed for interfaces are expected when large velocity contrasts exist and should be ignored. Globally, outside these hatched areas velocities can vary from ± 0.10 to ± 0.15 km s⁻¹ with rare excursions over ± 0.2 km s⁻¹.

4.4 MZ3 PSDM

A pre-stack depth migration and a residual moveout analysis were performed in order to convert MCS data from time to depth and verify the accuracy of the wide-angle velocity model. The Kirchhoff pre-stack depth migration procedure is described in detail by Schnürle et al. (2018), and the resulting section for MZ3 profile is shown in Fig. C1a. As for the pre-stack time migration (Fig. 3), the full sedimentary cover is well imaged up to the top of the SV1 layer. Below, the section is mostly transparent with only intermittent and tiny events within the SV1 layer. In the common-image gathers (CIGs: Fig. C1b), the top of Sv1 coincides with a strong low-wavelength reflector, followed by poorly coherent and under-migrated events. The semblance panels (Fig. C1c) are similarly well focused in the upper sediment and at the top of the SV1 interfaces, whereas it is blurred at greater depth and beneath the ALR.

5 Crustal structure of the Limpopo margin

The MZ3 final velocity model (Fig. 13) consists of 13 layers beneath the seafloor: 3 for upper sediments (S1–S3), 2 for the volcano-sedimentary sequence (SV1–SV2), 5 for the crystalline basement (G1–G5) and 3 for the mantle (M1–M3). In this section, we present the new constraints that the MZ3 MCS profile and velocity model bring to the crustal structure of the LM. We further include results from a tectono-stratigraphy analysis we made along key industrial MCS profiles in the area and from the other velocity models built within the scope of the PAMELA project along adjacent profiles. This allows us to give an overview of the crustal structure of the LM from top to bottom, including an interpretation of the stratigraphy of its upper sedimentary package, and to analyze the nature and deformation of its acoustic basement to finally focus on the crustal nature of the LM segments.

Figure 13 MZ3 final P-wave velocity model with region constrained by seismic rays highlighted. Blue lines indicates where wide-angle reflections constrain the model interfaces that are otherwise shown with black lines. Inverted black triangles mark OBS positions. Vertical exaggeration is 1:3.

5.1 Seismic stratigraphy of upper sediments

Our final velocity model (Fig. 13) shows two upper sedimentary layers over the NNV, increasing to three layers over the LM (Fig. 3). They are less than a kilometer thick combined around the ALR and progressively thicken as the seafloor and acoustic basement deepens toward the NE. This package shows globally low velocities between 2 and 3 km s⁻¹ top to bottom, reaching values slightly above 3 km s⁻¹ at the northeastern end of the line where it is the thickest (3.5 km). It is therefore largely made of uncompacted terrigenous sediments, except for close to the ALR where strong velocity contrasts and inversions suggest interleaved volcanic layers. Over the LM, sedimentary horizons are strongly deformed, showing typical features of a margin shelf and slope with incised valley or reworked sediments, such as mass transport deposits and contourites (Figs. 3, 14).

Figure 14 Compilation of line drawing of commercial MCS two-way travel time profiles superimposed by stratigraphic interpretation. The map in the bottom-left corner gives the locations of these profiles with respect to seismic profiles acquired during the PAMELA MOZ3-5 cruise. Also reported are drills locations that help to constrain ages of early sedimentary deposits. We also plotted velocity versus depth profiles at locations where profiles B and C cross the MZ3 velocity model.

On the southern MCP and offshore from Mozambique, several wells were drilled into this upper sequence and some of them were drilled up to the acoustic basement. Combined with commercial seismic profiles (Fig. 14), we performed a seismostratigraphy analysis to constrain the ages of the oldest deposits that span the syn-rift and early post-rift period that are of interest for this study. More precisely, we extended to the LM previous seismostratigraphy studies produced for the NNV (Baby et al., 2018; Schnürle et al., 2018). Profile A is the direct north–south prolongation along the LM of a profile presented by these authors. It connects to the Sunray-1 well (Salman and Abdula, 1995) located in the NW corner of the NNV (Fig. 14). We further analyzed two profiles (B and C) that are perpendicular to profile A and strike roughly E–W across the LM. Among them, profile B has an onshore portion that connects to Funhalouro-1 and Nhachengue-1 wells and Line SM-59 presented by Salman and Abdula (1995). Along all these lines, the horizons of late Jurassic to Neocomian formations, namely the Red Beds formation (Fm) and Maputo Fm (in green and yellow, respectively, in Fig. 14), can be confidently extended towards the LM from both the NNV (Sunray-1 well) and the MCP (Funhalouro-1 and Nhachengue-1 wells). As mentioned by Salman and Abdula (1995), these two formations unconformably cover the acoustic basement. Red Beds Fm is only observed filling local fossil depressions on the MCP and northern NNV (Fig. 14, green deposits). A possible similar deposit is interpreted at the eastern foot of the basement high on profiles B and C. On top of these restricted deposits or directly on the acoustic basement, the Maputo Fm appears widespread over both the NNV and the MCP. Most importantly, a continuous sequence of Maputo Fm (in yellow), which is also strongly eroded, can be drawn on profile B across the entire LM. It indeed passes clearly over to the east of the basement high where the first marine sediments (in blue in Fig. 14) are seen onlapping the formation. This constrains the age of the LM acoustic basement to pre-Neocomian time and before the onset of seafloor spreading on the MB.

Figure 15 Map of the main segmentation and boundaries within our study area resulting from the combined interpretation of PAMELA-MOZ35 deep seismic profiles. Pink lines are locations of line drawings of commercial MCS profiles shown in Fig. 14. The background shows the main geological units and structures on land: KC stands for Kaapval Craton, LB stands for Lebombo monocline, MSM stands for Mateke Sabi Monocline, MCP stands for Mozambique Coastal Plain, NNV stands for North Natal Valley, and ZC stands for Zimbabwe Craton. Offshore bathymetric contours in the LM are as follows: Limpopo margin, Mozambique Basin (MB) and South Natal Valley (SNV). The main geological features are as follows: Almirante Leite Ridge (ALR), Beira continental block or high (BH), Limpopo Fault (LF), Mozambique Ridge (MR), Mozambique Fracture Zone (MFZ), Tegula and Naude ridges, and Ariel Graben (AG). MOZ35 seismic acquisition in the NNV and LM is shown in red, with MZ3 profile highlighted with a bold red line.

5.2 Nature and deformation of the acoustic basement

At the base of upper sediments, the top of the acoustic basement, as identified on the MZ3 MCS line, consists of a set of high-amplitude parallel reflectors over the NNV and western LM changing to a rough surface to the NE end of the line (top of SV1 in Fig. 3). In the velocity model (Fig. 13), it is characterized by a large velocity jump (1–1.5 km s⁻¹) to reach velocities of 4.5–5 km s⁻¹. It is relatively flat over the NNV at about 2 km depth and deepens from 3.5–4 km in the LM to 6 km after 330 km model distance. Its perturbed topography is evidence of important deformations that are most probably controlled by highly dipping and near-vertical faults. On each side of the ALR they bound a few horsts and grabens, while over the LM they are highlighted by fan-like dipping reflectors and steep basement highs. The most prominent basement high, located at 310–340 km model distance, is clearly an isolated basement structure controlled on each side by steep faults. To its SW, the fault zone appears wider at the basement surface, suggesting a deep fault rising as a flower system (Fig. 3).

Internally, the acoustic basement shows a chaotic and incoherent signal in some places (60–120, 300–420 km, Fig. 3) but clear basement-parallel or dipping reflectors of variable amplitude in others (−60 to −20; 120–260 km). This mixed nature already suggests that the top of the acoustic basement is not the roof of the crystalline basement, at least up to the basement high at 330 km where it is more homogeneous. Its internal structure is, however, only partially resolved by the MCS. WAS data further constrain a continuous, 3 to 4 km thick unit made of two layers (SV1 and SV2 layers in green in Fig. 3b). They are, respectively, 1–1.5 and 1–2.5 km thick with velocities ranges of 4–5 and 5–6 km s⁻¹ on average, though slightly lower values are observed beneath the ALR (Fig. 13). From its seismic facies and velocities we therefore argue for the presence of a deep volcano-sedimentary unit that is continuous from the NNV towards the LM. This unit possibly ends at the basement high at 330 km model distance. Beyond 340 km model distance, our velocity model only constrains a single layer (a prolongation of SV1 for modeling purpose), which is 1.5 km thick with intermediate velocities of 5.0–5.5 km s⁻¹. From these velocities and its homogeneous, high-amplitude but chaotic facies, we instead suggest the presence of extruded volcanism only at the NE extremity of MZ3.

Commercial MCS data and companion PAMELA profiles complement and clarify the nature and extent of the deep volcano-sedimentary unit (SV1 and SV2). As for MZ3, all those profiles constrain a similar volcano-sedimentary unit (Moulin et al., 2020; Leprêtre et al., 2021; Schnürle et al., 2018) that ends oceanward in the form of a prominent basement high (Watremez et al., 2021). Along PAMELA MOZ3-5 WAS profiles (MZ1, MZ6, MZ7) to the southwest, an identical unit in term of thickness and velocity consisting of parallel and undulating horizons is seen striking through the entire NNV before it dips oceanward to the south or east (Moulin et al., 2020; Leprêtre et al., 2021; Schnürle et al., 2018). At the southern LM, along MZ3 crossing profiles MZ5 (Watremez et al., 2021) and C (brown color in Fig. 14), a broad basement high consists of internally gently dipping deformed but parallel reflectors. Further north, along MZ4 (Watremez et al., 2021) and profile B (Fig. 14), it is a set of landward-dipping reflectors (LDR) that flank the edge of the basement high. These LDR, however, connect with deep, wavy and parallel horizons seen extending from the eastern MCP. There is therefore a connection and a prolongation of a volcano-sedimentary unit from the MCP, through the western LM and up to the NNV that is further highlighted along profile A, where it is seen continuously all along the MCP coastline (Schnürle et al., 2018) (brown color in Fig. 14. All together, these observations reveal that an ubiquitous volcano-sedimentary unit exists over our study area from the southeastern MCP to the entire NNV (in green in Fig. 15). We therefore suggest that this unit forms a widespread pre-Neocomian basin that terminates to the east along the LM, where it is strongly deformed and associated with a prominent basement high and possibly also deeply rooted faults.

5.3 Crustal nature and segmentation

Below the volcano-sedimentary basin, the top of the crystalline crust is higher at 4–5 km depth in the NNV than in the LM where it reaches 6–7 km depth on MZ3 (Fig. 13). Considering the presence of a volcanic layer instead of a basin oceanward of the basement high, the top of the crystalline crust is thus deeper at 7–7.5 km depth. Along the entire MZ3 profile this interface marks a velocity jump to values above 6 km s⁻¹ that is responsible for a strong reflective phase observed on OBS records. Most importantly it is the place from which the velocity gradient of the crust becomes much lower around 0.044 km s⁻¹ km⁻¹ compared to the 0.25–0.3 km s⁻¹ km⁻¹ gradient within the upper acoustic basement.

Over the NNV (−50–110 km) the five layers that make the crystalline crust represent a 30 km thick unit with velocities between 6.0 and 7.3 km s⁻¹ top to bottom and a Moho lying at 34 km depth. Over the unthinned crust of the LM (up to 230 km) the same five layers represent a 24 km thick crust with a velocity range of 6.2–7.3 km s⁻¹ and a 31 km deep Moho. The strong reflectivity of crustal phases observed on OBS records over both areas is responsible for the numerous layers modeled within the crust and mantle and attest to their important magmatic content. Crustal thinning occurs over a distance of 100 km to the SW of the basement high where the Moho rises from 31 to 16 km deep (230–330 km). It coincides with a slight increase in velocity at the top of the crystalline basement (from 6.2 to 6.4 km s⁻¹) and a progressive decrease at its base (from 7.3 to 7.0 km s⁻¹). The velocity gradient thus remains stable and low (red curves in Fig. 16b), while its internal highly reflective character is still observed. Therefore, despite crustal thinning, the overall structure of the NNV and western LM are similar and consist of an upper volcano-sedimentary basin, a low-velocity gradient within the crystalline basement, strong magmatic intrusions and anomalously high velocities at the base of the crust.

These characteristics are also those observed and highlighted in velocity models produced by companion studies along other PAMELA Moz3-5 profiles that either cross the MZ3 profile at the LM (Watremez et al., 2021) or image the southeastern MCP and the NNV (Moulin et al., 2020; Schnürle et al., 2018; Leprêtre et al., 2021). Leprêtre et al. (2021) and Moulin et al. (2020) have made an extensive comparison of the crustal structure of the MCP and NNV with published velocity structures of thickened oceanic crust or continental crust from different tectonic settings. They conclude on the continental nature of the MCP and NNV crust, while the observed high velocities are explained by a its strong magmatic content. By comparing 1D velocity–depth (vz) profiles extracted every 10 km along MZ3 with vz profiles from MZ7 over both the MCP and NNV (Fig. 16a–b), we can observe the coherence in crustal structure between the MCP, the NNV and the LM despite the latter's progressive oceanward thinning compared to the normal continental crust (e.g., Christensen and Mooney, 1995). From the presence of a similar acoustic basement and similar crustal velocity profiles (vz), we argue that the MCP and NNV continental crust further extends to the east under the LM up to at least the location of the basement high (Fig. 15). This interpretation contradicts previous studies that either consider the entirety or parts of this region as floored by oceanic crust (e.g., Mueller and Jokat, 2019).

Figure 16 (a) The 1D velocity–depth (vz) profiles extracted along MZ3 final velocity model and color coded according to main domains and segments. (b) MZ3 vz profiles below the seafloor compared to MZ7 vz profiles of the MCP and NNV, and compiled vz profiles of continental crust for different settings (Christensen and Mooney, 1995). (c) Comparison of vz profiles below the acoustic basement of the corridor of anomalous crust extracted from MZ1 (Moulin et al., 2020; Schnürle et al., 2018), MZ3 (this study), MZ4 and MZ5 (Watremez et al., 2021) models with a compilation of vz from “normal” Atlantic oceanic crust (White et al., 1992). (d) Comparison of MZ3 vz profiles below the acoustic basement of the corridor of anomalous crust with vz for the oceanic crust of the MB extracted from MZ1 (Moulin et al., 2020; Schnürle et al., 2018), MZ4 and MZ5 (Watremez et al., 2021) and compared to the compilation of White et al. (1992).

On the northeastern side of the basement high, the crystalline crust of the LM is only 10 km thick in the MZ3 velocity model (Fig. 13). It is made of three layers with velocities between 6.5 to 7.0 km s⁻¹ top to bottom that show little internal reflectivity. While the two deepest layers present a low gradient (0.025 km s⁻¹ km⁻¹) similar to the adjacent continental crystalline basement (Fig. 16b), the upper layer has an intermediate velocity gradient (0.1 km s⁻¹ km⁻¹) between those of the acoustic basement above and the deep crystalline basement below (Fig. 16c). This anomalous trend is also observed on the “transitional crust” identified on MZ4 and MZ5 (Watremez et al., 2021) and east of the basement high on MZ1 velocity model (Moulin et al., 2020) as shown by the comparison of vz profiles from these profiles (Fig. 16c). Taken all together, they reveal a 60–80 km wide, N–S-trending corridor of anomalous crust isolated between the thinned MCP and NNV continental crust to the west and the MB oceanic crust to the east (see corridor of anomalous crust in Fig. 15).

The basement high therefore locates an important segmentation in the crustal structure of the LM (Fig. 15). Not only does it coincide with the eastward termination of the volcano-sedimentary basin and a possible change in nature of the acoustic basement but it also marks a profound modification of the crystalline crust. We can further notice that uppermost mantle velocities appear normal and stable at 8 km s⁻¹ on each side of the basement high, while a slight velocity decrease to 7.9 km s⁻¹ (Fig. 13), which is more pronounced in the MZ4 and MZ5 velocity models (Watremez et al., 2021), is present beneath the feature. This might suggest that the basement high is the surface expression of a vertical frontier that is deeply rooted in the lithosphere. As such, it is comparable to the MFZ that bounds the opposite side of the corridor and separates it from the oceanic crust of the MB. We therefore propose to refer to this zone of strongly localized deformation as the Limpopo Fault (LF, Fig. 15).

To summarize, the LM is a N–S elongated margin cut by two major fault zones that segment its crustal structure into a western continental part, a central corridor and an eastern oceanic part (LF, Fig. 15). The continental domain is the direct eastward prolongation of the continental crust that floors the MCP and NNV and is capped by a pre-Neocomian volcano-sedimentary basin. Defining the nature of the crust within the corridor remains a challenge. Vz profiles (Fig. 16c–d) reveal a possible mixed composition. On one hand, a 3–4 km thick upper part with a typical signature of oceanic crust is caused by the presence of extruded volcanics. On the other hand, a lower part that still preserves the low gradient characteristics of adjacent continental crust is also present. Its crustal thickness is also extremely variable (5 to 15 km) and appears thicker to the south (MZ1) than to the north (MZ3, MZ4, MZ5: Fig. 16c). Vz profiles of the MB oceanic crust located just east of the MFZ also shows some changes, with a velocity signature close to the corridor to the north (where it is the oldest) and a signature that shifts toward a typical oceanic crust to the south (Fig. 16d).

6 Discussion

6.1 Geodynamic implications

Deep crustal velocity models produced from the PAMELA-MOZ35 seismic experiment have revealed the presence of a 40 and 30 km thick crust beneath the southern MCP and offshore NNV, respectively (Moulin et al., 2020; Schnürle et al., 2018; Leprêtre et al., 2021). They also show a continuous and mildly deformed volcano-sedimentary cover over these two entities. These findings are strong evidence in favor of the continental nature of the MCP/NNV as suggested previously by seismological and magnetic studies (Domingues et al., 2016; Hanyu et al., 2017). This work further details the crustal architecture of the area, delimiting the eastern edge of this continental unit along the LM. The presence of an homogeneous basin over a thick continental crust attests to an intra-continental depositional process (Moulin et al., 2020) rather than a rifted margin process (e.g.,  Cox, 1992; Klausen, 2009; Watkeys, 2002). Our tectono-stratigraphic analysis only places an upper bound to the formation of the basin in pre-Neocomian. So far it has not been dated more precisely despite the fact that its upper volcanic layer was reached by several wells (Flores, 1973). We can thus speculate that it formed during the Karoo phase, which affected the entire African continent in late Palaeozoic–early Triassic (Daly et al., 1991). However, this would leave large, 30 Myr sedimentary gaps between extrusions of Karoo volcanics, which might compose the roof of the basin, and the deposition of the late Jurassic to Neocomian Red Beds formation and Cretaceous Maputo formation. Therefore, another hypothesis is that it formed in the late Jurassic either before or contemporaneously with the formation of the adjacent MB as suggested by (Li et al., 2021). The high magmatic content of the basin, as well as evidence of strongly intruded continental crust, suggests that magmatism may have overloaded the crust and created the necessary vertical subsidence (Moulin et al., 2020; Tozer et al., 2017).

Our study also emphasizes the eastward prolongation and termination of the volcano-sedimentary basin along the LM, where it is strongly deformed by deeply rooted vertical faults. At depth, it is also the place where important crustal thinning and segmentation is evidenced (Fig. 15). Therefore, the LM appears naturally as the place where rifting is localized to accommodate the opening of the MB rather than over then entire MCP and NNV. In previous scenarios based on “tight” fit kinematic framework in which Antarctica partly overlaps Africa, rifting was either postulated to concentrate along the Lebombo monocline or beneath the MCP and NNV depending on whether the area was interpreted as oceanic or as a volcanic rifted margin (Cox, 1992; Klausen, 2009; Leinweber and Jokat, 2012; Mueller and Jokat, 2019; Watkeys, 2002). Neither of these hypothesis are supported by these new observations.

Overall, this means that the MCP/NNV must be excluded from the Africa–Antarctica corridor (AAC in Leinweber and Jokat, 2012; Mueller and Jokat, 2019), and a “looser” plate fit must be adopted in East Gondwana kinematic reconstruction (Moulin et al., 2020; Thompson et al., 2019). Such a framework excludes any initial rifting phase with normal (E–W) or oblique (NW–SE) plate movement (e.g., Cox, 1992; Reeves et al., 2016), which in our opinion has never been clearly evidenced or described. Invoking the orientation of magmatic dike swarms to attest stress field direction (Mueller and Jokat, 2019; Reeves et al., 2016) is highly speculative as they may be strongly controlled by inherited lithospheric discontinuities (Jourdan et al., 2006). Similarly, the presence of a wide crustal necking zone cannot solely justify normal or oblique rifting (e.g., Vormann et al., 2020). Indeed the LM itself shows such characteristics, but given our preferred and alternative geodynamic framework it was affected instead by strike-slip or slightly transtensional rifting following a continuous N–S direction of plate motion during the opening of the MB.

6.2 Strike-slip rifting along the Limpopo margin

Along the LM, strike-slip rifting is emphasized by the Limpopo Fault (Fig. 15). It forms a zone of deeply rooted and strongly localized deformation rising to the surface as a flower structure (Figs. 3, 14). Basement uplift is evidenced by the presence of a prominent high along the fault with deposition of pre-Neocomian Red Beds formation on either side and signs of erosion of the Maputo Sands formation that covers it. Taken together, these indicators reflect strike-slip strain along the LF that lasted up to the early Cretaceous. Despite a thick sedimentary cover and strong magmatic content that both affect the definition of seismic profiles, we identified a long wavelength wavy deformation within the volcano-sedimentary basement that looks like folding, as well as a few horsts and grabens (Figs. 3, 14). These are located at the western edge of the LM and might reflect some partitioning in the deformation of the upper crust with a slight normal component. Uplifts, en echelon pull-apart structures and shear folding are all typical features that have so far been observed along highly oblique or strike-slip margins (e.g., Mascle and Blarez, 1987), but those contexts are more generally favorable to strain partitioning and a wide range of deformation features (Brune, 2014; Teyssier et al., 1995).

Besides upper-crustal brittle deformations of the LM during its rifting, our observations point toward important crustal reworking at depth accompanied by volcanic extrusions. From the thinned continental crust structure to the west of LM and the lower continental crust velocity signature of the anomalous corridor (Fig. 13 and the MZ1 model in Moulin et al., 2020; Schnürle et al., 2018), we infer that oceanward flow of the MCP and NNV lower crust may have “fed” the corridor. The corridor's crystalline basement is 10–12 km thick in MZ3 (Fig. 16c) in the northernmost LM, varies from 8 to 14 km thick offshore from the MCP (MZ4 and MZ5 in Fig. 16c) and thickens to 15–20 km east of the NNV (MZ1 in Fig. 16c). Such variation in thickness may be evidence of the “boudinage” of the lower crust originating from its ductile shearing and flow (Clerc et al., 2018; Gernigon et al., 2014; Loureiro et al., 2018). Strong volcanism has accompanied this process, as suggested by the “mixed” composition of the corridor and time constraints given by the deposition on top of these volcanics of the Maputo Sands that preceded the first deep marine sedimentation (Fig. 14).

Figure 17 Evolution of the Limpopo margin (LM). (a) Kinematic of the Antarctica plate with respect to the African plate fixed according to (Thompson et al., 2019). (b, c, d) Zooms on the evolution of the LM with major features annotated before rifting (b), at breakup (c) and after the opening of the Mozambique Basin along the LM (d). Magnetic anomalies are from Mueller and Jokat (2019). LB stands for Lebombo monocline, MSM stands for Mateke Sabi Monocline, MCP stands for Mozambique Coastal Plain, NNV stands for North Natal Valley, ME stands for Maurice Ewing bank, BH stands for Beira high, OZD stands for Offshore Zambezi Depression, AGA stands for Antarctic Gravity Anomaly, GC stands for Grunehogna Craton, LF stands for Limpopo Fault, MFZ stands for Mozambique Fracture Zone, AG stands for Ariel Graben and TJ stands for Triple Junction.

Lower continental crust flow and shearing, eventually leading to exhumation (Aslanian et al., 2009), has been inferred in many places including at the offshore eastern Canada and Arctic margins (Gernigon et al., 2014; Reid and Jackson, 1997), along the South Atlantic Brazilian (Evain et al., 2015; Loureiro et al., 2018) and western African margins (Clerc et al., 2018; Moulin et al., 2005), offshore of Mozambique (Senkans et al., 2019), and in the Mediterranean's Gulf of Lions and Sardinian conjugate margins (Afilhado et al., 2015; Moulin et al., 2015; Jolivet et al., 2015). It appears as a common process of rifted continental margin independently of its tectonic context. In fact it has been shown numerically that lower crust may flow across both divergent margin segments (Huismans and Beaumont, 2011) and strike-slip segments (Le Pourhiet et al., 2017; Reid, 1989). Usually it can be seen that the weak lower crust may favor such an exhumation process. This was undoubtedly the case for the southern Mozambique lithosphere, which on the verge of its dislocation was affected by a large thermal anomaly (possibly Karoo) responsible for its anomalous magmatic content. The same anomaly might have favored coeval partial melting and extrusion of volcanic material to form the volcanic basement of the corridor.

6.3 Rift to drift evolution and vertical movements

In the alternative East Gondwana “fit” proposed by Thompson et al. (2019), the continental domain to the southwest of the Astrid Ridge on the Antarctica plate is the conjugate of the African corner made by the MCP and the Beira High (Fig. 17a–b). Further south, the Grunneghona craton is facing the NNV. The exact seaward limit of the Antarctica continental crust is inferred at the location of a strong positive free-air gravity anomaly (Mueller and Jokat, 2019; Scheinert et al., 2016) since transitional crust is observed further offshore in seismic data (Jokat et al., 2004). At this stage, the Antarctica plate overlaps the corridor, which develops progressively during rifting (Fig. 17b) and before oceanic spreading in the MB starts from Chron M25 (155 Ma) onwards (Fig. 17c–d). As discussed above, we infer an intra-continental sag phase for the development of the volcano-sedimentary basin over the MCP and NNV. Deep crustal magmatic intrusions may have led to the subsidence of the area without major plate movements (Moulin et al., 2020). On the other hand, rifting along the LM proceeds from an almost N–S plate motion responsible for the opening of the MB (Fig. 17b–c). The LM develops as a wide shear zone along the eastern margin of the MCP and NNV. The LF acts at this time as a major strike-slip fault while easing the decoupling between the upper and lower crust at depth. In the upper crust we infer that strain partitioning is responsible for shear folds, small pull-apart basins and grabens within the volcano-sedimentary basement. Deeper extension is accommodated by lithospheric flow leading to lower crustal thinning of the eastern fringe of MCP and NNV continental crust and its oceanward exhumation. This was accompanied by the extrusion of a large amount of volcanics over the corridor. Further west, deep seismic acquisition revealed a transitional domain capped by a thick volcanic layer (the Explora Wedge, e.g., Jokat et al., 2004) off the Antarctic margin, which might attest to an equivalent process, albeit along a divergent segment. On the African side, extension was initially focused in the offshore Zambezi depression before migrating to the south of the Beira High, isolating this continental block (Mahanjane, 2012). The nature of the crust to the north of the Beira High is also deemed to be transitional with syn-rift magmatism reported (Mueller et al., 2016; Mueller and Jokat, 2019).

According to our stratigraphic analysis, the scarce Red Beds formation continental sedimentation was followed by the widespread deposition of the shallow marine Maputo Sands formation over both the MCP and NNV and the LM. Only then was the eastern extremity of the corridor covered by the first deep marine sediments. Lithospheric and/or asthenospheric flow may have sustained the area to a relatively high level during continental rifting (Reid, 1989) first allowing the continental Red Beds formation deposits, followed by a slight subsidence responsible for the first shallow marine incursion. Minor vertical movements must have occurred along the LF only due to strike-slip shearing. It is some time after the breakup that larger differential vertical movement took place on either side of the LF between the continental domain and the corridor, accentuating the basement high and limiting the westward incursion of the first deep marine sedimentary horizons.

As in previous studies (see  Basile, 2015, and references therein), we recognize that a stage of continent–ocean interaction followed continental rifting and its break-up. Analyzing the possible interplay between the corridor and the mid-oceanic ridge requires, however, a careful analysis of later stratigraphic horizons across and all along the margin that is out of the scope of our study. Figure 17c to d illustrates the drift of the mid-ocean ridge axis through time based on identified magnetic chrons (Mueller and Jokat, 2019, and references therein). Valanginian (Fig. 17d) corresponds to a period of major kinematic reorganization in the area. While the Patagonian block started its southwestward drift, a triple junction initiated that is responsible for enhanced magmatism and the formation of the Mozambique ridge (Fischer et al., 2017). At this time the mid-oceanic ridge axis reached the position of the Ariel Graben. Active deformation thus ceased along the LM, which entered into a passive stage and progressively acquired it present morphology.

7 Conclusions

This study focuses on the MZ3 combined MCS and WAS profile acquired across the LM within the scope of the PAMELA MOZ3-5 project. A P-wave velocity model of the deep crustal architecture of the margin is produced and interpreted with additional inputs from a tectonostratigraphic analysis of industrial seismic profiles and adjacent velocity models from companion studies (Moulin et al., 2020; Schnürle et al., 2018; Leprêtre et al., 2021; Watremez et al., 2021). We evidence the eastern prolongation, thinning and termination of the MCP and NNV continental crust under the LM. We further highlight its segmentation emphasized by the presence of a N–S-oriented corridor of anomalous crust bounded by deeply rooted faults, i.e., the Limpopo fault to the west and the MFZ to the east. Consistent with recent plate kinematic reconstruction of East Gondwana breakup, we infer that strike-slip or highly oblique rifting occurred along the LM. This led to the thinning of the continental crust beneath the eastern edge of the MCP and NNV, while the observed peculiar “mixed” crustal composition of the corridor is explained by lower crust flow with additional magmatic inputs. Finally, we show that only moderate subsidence occurred during the intra-continental rifting and early continent–ocean interaction stages. It is only after the break-up that decoupling along the LF led to differential vertical motion on each side of the fault with uplift of the continental domain while subsidence affected the corridor.

Appendix A

Table A1 Refracted (left) and reflected (right) phase names for MZ3 profile, number of explained events (npts), residual mean squares (Trms, in seconds), and normalized chi-squared values.

Table A2 OBS names and distance along the MZ3 model, phase propagation direction, number of explained events (npts), residual mean square (Trms, in seconds), and normalized chi-squared value.

Appendix B: MZ3 uncertainties estimation

Figure B1 MZ3 final velocity model with parameter selection (white hexagons) for VMONTECARLO. The top panel shows depth parameters, the middle panel shows the top of the layer velocity parameters and the bottom panel shows the base of the layer velocity parameters. Interfaces are indicated by black lines.

Appendix C: MZ3 PSDM

Figure C1 Pre-stack depth migration of MCS MZ3 profile. The left column shows the SW portion, and the right column shows the NE portion. (a) PSDM record section. The intersections with the MOZ35 dataset are indicated by a red line. Vertical exaggeration is 1:10. (b) Residual moveout. (c) Semblance plot of the RMO. Common image gathers are spaced every 7.5 km. Vertical exaggeration is 1:15. Model interfaces are represented with continuous lines.