Crustal Structure beneath the Precambrian Cratons of Gondwanaland and Its Evolution Using Teleseismic Receiver Function

Cratons are stable portions of the Earth’s crust that have undergone very little deformation for prolonged time periods [1]. The Precambrian cratons have survived the cycles of merging and rifting of continents over geologic time, making them unique. In order to understand how the cratons formed and evolved, investigation of the lithospheric structure is essential. During the Gondwana period, the landmasses of South Africa, Antarctica, Australia, Madagascar, Sri Lanka, and the Indian subcontinent were all joined near the South Pole, making the Gondwana supercontinent an accretion of several cratons (Figure 1). The breakup of Gondwanaland during the Mesozoic led to the gradual rifting of these different cratons over geologic time at different rates. Sri Lanka had also moved along with the Dharwar craton and shared similar crustal properties with its Gondwana neighbours [2]. Cratonic roots are supposed to have a very strong thick lithosphere with thick crust [3]. But after the Gondwana breakup, the cratons have migrated along with the continents at different speeds, different directions, and some travelling very large distances. One of the important concerns is to what extent the crustal properties were affected during the continental evolution history after the Gondwana breakup. Another concern is whether there is a link between evolution of cratons and their seismic structure.

Models of crustal growth can be broadly categorized as crust based and mantle based. The crust-based models focus on cumulative crustal material at the present day, but the mantle-based models also shed light on crust that has been lost to mantle by reworking, subduction, or delamination, hence focusing on net crustal growth [4]. These studies are also aimed at estimating the formation age and how crustal processes lead to deviations in surface age, formation age, and average age of the crust. Cratons are the oldest regions on the earth preserving old crustal material, and because of this fact, investigating the crustal architecture and crustal properties with respect to the age of the cratons is crucial in understanding crustal evolution. Opinion remains divided in the scientific community with respect to crustal properties and ages of crustal material. One of the models by Smithson et al. [5] indicates that the crustal thickness and mean seismic velocity are proportional to age (except young orogens), implying greater crustal thickness should be observed in Archean crust compared to Proterozoic aged crust. An alternate view by Durrheim and Mooney [6] suggested that the velocity at the base of the Proterozoic crust is high (⁠$Vp>7.0$⁠) and is thicker (40-50 km), when compared to the Archean crust (27-40 km). Nelson [7] suggested that repeated episodes of the basaltic underplating cause crustal thickening. Rudnick and Fountain [8] found similar crustal thickness of Archean and Proterozoic. Arndt and Goldstein [9] had suggested that transportation processes between the lower continental crust and the upper mantle are potential contributors in building of the continental crust. It is still uncertain how crustal formation, growth, and composition have varied during the Precambrian; these open questions are addressed in this study.

The Moho or the Mohorovicic discontinuity defines the base of the crust and is the outermost seismic discontinuity in the Earth [10]. It is generally marked by an increase in velocity of P wave from <7.6 to >8 km/s [11]. Petrologically, it is called the crust-mantle transition zone, where the composition changes from mafic rock to ultramafic rock. Even though both terms are used simultaneously, studies suggest that they might vary [12]. According to Griffin and O’reilly [13], thick mafic roots may convert to the high-density mafic rock eclogite at base of crust. Eclogite has similar seismic velocities to the mantle, and hence, the crust-mantle boundary may not be observable as an abrupt increase in seismic velocity. Therefore, it is crucial to image the crustal properties and study the Moho nature. If, in a region, the lower crust is felsic and is underlain by a mafic upper mantle, then the nature of the Moho will be seismically sharp. On the other hand, if a mafic layer with high velocity is present at the base of the lower crust (as can be caused by intrusion from basaltic underplating) and is underlain by mafic upper mantle, then the Moho is seismically diffusive or gradational in nature.

Thus, we see that there is a link between evolution and seismic structure. Study of crustal structure helps in narrowing down on the competing hypothesis for crustal evolution. Repeated episodes of basaltic underplating have the possibility to manifest as a high velocity mafic layer at the base of the crust and could lead to crustal thickening. Study of crustal structure of different ages is also one of the best methods to test the different models of crustal evolution. The best areas to test this would be old crustal material that has not been deformed for a prolonged period of time and have survived the cycles of merging and continent rifting, the Precambrian cratons.

In this study, crustal properties have been imaged, comprising crustal thickness, crustal composition (⁠$Vp/Vs$ ratio), and nature of the Moho of the Precambrian cratons of Gondwanaland, using receiver function modeling. The technique of a single method (receiver function) used in this study for determination of crustal structure beneath the different areas of Precambrian cratons of Gondwanaland contributes to the subject of crustal evolution. It is worth investigating the crustal properties across the Precambrian of different ages and investigates the link between seismic structure and evolution. In order to obtain the detailed shear velocity structure along with $Vp/Vs$ ratio variation with depth beneath the stations in different cratons, inversion of receiver function is needed. In this study, we obtain the shear velocity structure (⁠$Vs$ and $Vp/Vs$ with depth) from shallow to the Moho depth using two different methods, the $H−κ$ stacking and neighbourhood algorithm. Some of the previous studies in cratons just used the $H−κ$ stacking approach, which is dependent on a priori P wave velocity (⁠$Vp$⁠) information, which may not be correct every time. Wrong selection of $Vp$ may give erroneous $H$ and $κ$ values. Also, $H−κ$ stacking provides crustal thickness and average crustal $Vp/Vs$ ratio. Variation of $Vs$ and $Vp/Vs$ has important implications in terms of knowing any intracrustal discontinuities, or deep rocks in the crust, presence of fluids, etc.

Study area consists of Dharwar craton in India; Pilbara, Yilgarn, and Gawler cratons in Australia; Tanzania, Kaapvaal, and Zimbabwe cratons in Africa; Antananarivo craton in Madagascar; and Napier complex in Antarctica. The brief geology of all the cratons are discussed below.

It is located in the southern part of India (Figure 2(a)). The crust of the Dharwar craton dates back to 3.4 Ga. The different lithological units are the Eastern Dharwar craton (EDC) and the Western Dharwar craton (WDC) that are separated by the Chitradurga schist belt (CSB) [14]. The Kaladagi basin of Proterozoic age borders the WDC in the north. Granulites and low-grade gneiss are present in the South of WDC, and greenstone is present in the North. To the west, the Western Ghats surrounds the WDC. There is a signature of thrusted contact between EDC and the Proterozoic Cuddapah Basin. On the eastern side, the Eastern Ghats surround the EDC. The Dharwar craton slowly transitions into the Neoarchean metamorphic terrane of the Southern Granulite Terrane in the south.

The craton lies in Western Australia (Figure 2(b)). The Proterozoic Capricorn Orogen in the south separates it from the Yilgarn craton. Zircon crystals having an age of 3.8 Ga are exposed in the early Archean craton. Lithologically, it is divided into Western Pilbara terrane, Eastern Pilbara terrane, Kurrana terrane, and the Hamerseley Basin [15]. The Elizabeth Hill belt was created when the eastern and western part of the Pilbara craton collided. Granitic complexes are present in the southeastern part of the Kurrana terrane, and granite greenstone units are present in the western and eastern part. The Neoarchean Hamersley basin is present in the south and consists of volcanic sedimentary rocks. Tonalite-trondhjemite-granodiorite (TTG) rocks, greenstone belts, and dome and keel structures are also present in the Pilbara region that makes it unique. It has been suggested that vertical tectonics operated in the Pilbara, and due to partial convective overturn, the dome and keel structures were formed [16].

This craton (Figure 2(b)) and the Pilbara craton are separated by the Capricorn orogen, which is formed during Paleoproterozoic collision between the two cratons. The age of the oldest zircons dates back to 4.2 Ga [17] during the Hadean period. It consists of many terranes—Murchison province, Western gneiss terrane, Eastern Goldfield province, and the Southern cross province [18–20].

The Gawler craton is situated in the southern part of Australia (Figure 2(b)). It comprises dominantly volcanic and sedimentary rocks. The age of the oldest granites in the craton dates back to 3.0 Ga. Tectonic activities in the Late Archean and early Proterozoic have led to the development of the craton. Felsic magmatism was dominant during the late Archean, and rifting was dominant initially during Proterozoic.

Tanzania craton is located in the Eastern part of Africa (Figure 2(c)). The lithology consists of different units made of granites in the north and gneiss and granites in the South separated by the Dodoma Greenstone belt. The Mozambique belt runs along the eastern side of the craton and the Usagaran belt lies to the southeast. The Mozambique belt was created during Pan-African orogeny as a suture zone during the formation of Gondwanaland. The Tanzanian craton, which is formed due to the overriding plate of the subduction zone, created the Usagaran belt [21, 22]. Greenstone belts, TTG (Tonalite-trondhjemite granodiorite rocks), and gneisses are exposed on the craton [23]. Xenolith and geophysical studies indicate that a rigid and deep lithospheric keel is present beneath the craton [24–26]. Presence of a superplume beneath the craton has led to the dynamic topography of the craton [27].

It is located in the Southern part of Africa (Figure 2(d)). Similar rocks from this craton and the Pilbara suggest that they were part of the Vaalbara supercontinent [28]. Detrital zircons from the Gneiss complex in Swaziland date back to 3.7 Ga [29]. The lithology consists of granites, greenstone belts, and gneisses. The Bushfield igneous province is present in the North, along with the Limpopo belt that joined it with the Zimbabwe craton during late Archean metamorphism. Komatiites are present in the Barberton greenstone belt in the east.

The Limpopo orogenic belt towards the Southern part of this craton separates it from the Kaapvaal craton (Figure 2(d)). It also consists of granites, greenstone belts, and gneisses. Lithologically, the Zimbabwe craton joined with the Kaapvaal craton in the Late Archean, and they are together called the Kalahari craton. The combination of sutures of Tokwe segment to the south and the Rhodesdale fragment in the north formed this cratonic block. Smith et al. [30] has found 3.7 Ga gneiss in the Tokwe segment. It has been suggested that the craton formed by arc accretion method [31]. In the northwest, it is surrounded by the Magondi belt and the Zimbabwe belt borders to the north of this craton.

It is located in central Madagascar (Figure 2(e)). The oldest rocks in the craton consist of migmatite gneiss that dates back to the Neoarchean age. The northern Sofia group and the southern Vondrozo group make up the supracrustal units of the craton. Intrusive origin migmatite gneisses are also present [32]. Three large synformal belts consisting of mafic schist and gneiss constitute the Tsaratanana complex in this craton. The rocks of the craton are of high metamorphic grade and complex in structure, and because it was formed in Neoarchean, it was subjected to magmatism, deformation, and metamorphism during the Proterozoic [33].

The Napier complex lies in Enderby Land in Eastern Antarctica (Figure 2(f)). Before the breakup of Gondwanaland, the Dharwar craton of India shared a common boundary with this craton [34]. Various tectonothermal events, high-grade metamorphism, and plate tectonics led to Napier orogeny and formed the cratonic nucleus. Subsequent years of erosion and tectonic deformation have exposed the lower metamorphic core of the ancient mountains. The oldest tonalite and granitic gneiss record an age of 3.8 Ga [35]. Younger granulite facies rocks derived from younger granites and granodiorites and deformed gneisses from volcanic and sedimentary protoliths are also present [36]. In the late Archean, ultrahigh-temperature metamorphic conditions have led to their deformation. High-temperature metamorphism in the early Archean at 24-40 km crustal depth at 1050-1120°C was followed by sustained cooling periods at 20-30 km depth [37]. A subglacial depression separates the remaining Antarctica from this block.

The teleseismic events (epicentral distance in between 30° and 95°) for the different cratons were chosen such that all have a magnitude greater than or equal to 6. All the data was acquired from the Data Management Centre (DMC) of the Incorporated Research Institutions for Seismology (IRIS). A total of 52 seismic stations are used in our analysis for different cratons. In our $H$ and $Vp/Vs$ ratio calculation, we have used all freely available data from all the stations of Gondwanaland. If for a craton, two or more stations have similar $H$ and $κ$ estimates, we have taken representative station(s) to estimate the crustal parameters. The period of earthquake data in Dharwar was from 2006 to 2007 using instrument Streckeisen STS1 with sample rate of 20 samples/s. The data was from 2015 to 2017 with a sample rate of 40 samples/s in the Australian cratons using the instrument KS2000 in Pilbara, Guralp CMG3T in Yilgarn, and Guralp CMG3ESP in Gawler. In Tanzania craton, the earthquake data was extracted for the year 1994 with a rate of 20 samples/s that used instrument STS2 in stations MTOR, SING, and PUGE and with instrument cmg3esp in station MBWE. In the GETA station of Tanzania, the earthquake data was from 2011 to 2015 with a rate of 40 samples/s using the instrument trillium compact. In the DODT station of Tanzania, the data was from 2009 to 2015 with a rate of 40 samples/s using the instrument Streckeisen STS-2. In Kaapvaal craton, the period of earthquake data was from 1997 to 1999 using the instrument STS2 with a rate of 20 samples/s. In Zimbabwe craton, the period was from 1997 to 1998 in stations SA72, SA74, SA75, and SA77 and from 1997 to 1999 in station SA76, all using instrument STS2 with a sample rate of 20 samples/s. In Antananarivo craton, the period of earthquake data was from 2007 to 2010 in station ABPO with instrument Streckeisen with a rate of 20 samples/s and from 2012 to 2013 in stations BARY, BATG, BITY, and VINA using instrument Streckeisen with a rate of 40 samples/s. In Napier complex, the earthquake data was from 2013 to 2016 using the instrument Streckeisen with 40 samples/s. The coordinates of the stations in the cratons along with the networks are given in Table 1.

The internal structure of the Earth near the receiver was estimated by modeling the time domain receiver function, which is calculated by deconvolving the horizontal component with the vertical component [38–40]. When plane waves interact with a discontinuity (say, the Moho), a part of P wave gets converted to an S wave (denoted as Ps), and other parts undergo multiple reflections (multiples) in between the Earth’s surface and that discontinuity (denoted by PpPs and PpSs). To compute the receiver function, ,first, the TauP-Toolkit developed by Crotwell et al. [41] has been used to mark the P arrivals in the waveforms with reference to the IASP91 velocity model. The mean and trend are removed from the waveforms, and 150 s of the waveforms (from 30 s before P arrival time to 120 s after P arrival time) is selected for further processing. The radial and transverse components are then obtained by rotating the north-south and east-west components with respect to a great circular arc. The method of iterative time domain deconvolution [42] has been used to obtain radial and transverse receiver function by deconvolving the radial and transverse components, respectively, with the vertical component. Lowpass Gaussian filter of width (⁠$Gw$⁠) 2.0 (corner frequency ~1 Hz) was applied to the receiver functions. The receiver functions with more than 80% fitting are used in this study. The Ps arrival times are increased, and their multiples are decreased as ray parameters increase. So, move-out correction of the receiver functions has been done for individual stations using the IASP91 velocity model and reference epicentral distance of 67°. Each station consists of many good quality receiver functions; the minimum number is kept to 5.

Figure 3 shows the receiver functions of one representative station from each of the Archean cratons. Variations of Ps, PpPs, and PpSs are also shown in the figure with different colors, which can provide initial information of crustal thickness beneath each station. The Ps arrival time for the station HYB of Eastern Dharwar craton, India (Figure 3(a)), and PSA00 of Pilbara, Australia (Figure 3(b)), is observed at ~3.6 s; however, their multiples (PpPs and PpSs) have different arrivals. For the other two cratonic blocks of Australia, Yilgarn (PSA00) and Gawler craton (MULG), the Ps arrival times are observed at ~4.5 s and ~5.5 s, respectively. This suggests a thicker crust beneath the Gawler craton, compared to Pilbara and Yilgarn. Pilbara may have the thinnest crust in Australia. Based on the Ps arrivals of receiver functions estimated for DODT of Tanzania (Figure 3(e)), SA38 of Kaapvaal, and SA74 of Zimbabwe, which is ~4.5 s, with PpPs arrival ~15-16 s and PpSs arrival ~19-21 s; we can say crustal thickness may be similar for these cratons. In Antananarivo craton (ABPO) and Napier complex (MAW), the Ps, PpPs, and PpSs arrival times are observed to be ~4.5 s, ~17 s, and~22-23 s, respectively, suggesting similar crustal thickness in both the cratons.

The average P wave velocity of the crust is considered 6.4-6.5 km/s for the calculation of $t1$⁠, $t2$⁠, and $t3$⁠. Vp was varied between 6.1 and 6.6 km/s to calculate the uncertainties in $H$ and $κ$ for the error analysis, along with changing the combination of weighting factors. The uncertainties varied from 3.0 km for $H$ and 0.5 for $Vp/Vs$⁠. The $H−κ$ result of station PSAD3 in Pilbara craton is shown in Figure 4.

The $H−κ$ method calculates the crustal properties by assuming an average $Vp$⁠, which may not be the same for all the stations. Hence, only calculation of Moho depth and an average $Vp/Vs$ ratio is possible. To find out $Vs$ and $Vp/Vs$ ratio variation with depths, inversion modeling is used. Some other previous studies used tomography or gravity inversion for large-scale regions; hence, the resolution for the crust could be low. Inversion of receiver function using the neighbourhood algorithm method is a nonlinear method free from a priori information, which makes the estimations more robust in the study. To reduce the effect of nonuniqueness and nonlinearity, which are common problems in receiver function inversion (Ammon et al., 1990), Sambridge [44] developed a fully nonlinear derivative free direct search approach, the neighbourhood algorithm (NA) for finding models of acceptable data fit. Through the modeling, we can get an ensemble of models that preferentially sample the good data-fitting regions of parameter space, rather than seeking a single optimal model. Linearized inversion is strongly dependent on the initial model. However, NA is less sensitive to the initial model, and in most cases, it converges to the global minimum. We also observed that for different initial models, NA produces a similar final model by only changing either the number of iterations and/or number of samples/resamples.

The final model is selected from the set of 1000 best models, based on the minimum misfit between the synthetic and observed receiver functions. Misfit is the measure of the discrepancy between the true and predicted RF from an Earth model. Minimum misfit does not guarantee the best model, and in many cases, it provides an unrealistic Earth model. We have given preference to the model providing minimum misfit (may not be the least) with a realistic Earth model. Detailed velocity modeling at depth is needed in order to explore the crustal nature. In this method, the lithosphere consists of 6 horizontal layers, namely, sediment, basement, upper, middle, and lower crusts, and upper mantle, with 4 model parameters, the $Vp/Vs$ ratio, the thickness of layer, and the shear wave velocities at the uppermost and lowermost of the layer. This makes it a 24- (⁠$6×4$⁠) dimensional parameter space. In this method, synthetic receiver functions are constructed for $ns$ models that were created for every iteration. An objective function is defined based on the $L2$ norm of the difference between observed and synthetic receiver function, which is considered the misfit. Then $nr$ models were selected based on lower misfit values. Then again, $ns$ new models are created based on random walk through the Voronoi cell of $nr$ models. In this way, new models are produced in areas where the fitting is good. The fitting and inversion result of station PSAD3 in Pilbara craton, Australia, is shown in Figures 5(a) and 5(b), respectively. 40,000 velocity models that have been generated by iterating the neighbourhood algorithm 400 times are shown in Figure 5(b), and the original stacked receiver function is compared with the synthetic receiver functions corresponding to each model. Best fitted (colored section) 1000 models are selected based on the minimum misfit between the synthetic and observed receiver functions. The red line represents the best fitted model. The observed stacked receiver function (black line) along with the best fitted synthetic receiver function (red line) bound by $±1σ$ error bounds (black dashed line) is shown in Figure 6 from one station in each craton. The Moho could be observed as a sharp step jump in the shear wave velocity inversion which is easy to identify, or as a gradational change in velocity. In the latter case, it is picked up where the velocity crosses 4.2-4.3 km/s (Figure 7). The same process is followed for all other stations (included in supplementary available (here)). The uncertainties were calculated from the best fitting models. The uncertainties for $κ$ and $H$ are shown in Figure 8.

Table 1 incorporates the $H−κ$ and NA results from the stations used in this study. The $H−κ$ results, fit of observed and synthetic receiver functions, and NA results from one station in each craton are shown in Figures 6, 7, and 9, respectively, and the remaining results from the other stations have been illustrated in supplementary and are documented in Table 1. The crustal thickness and $Vp/Vs$ ratio estimated using both the methods (⁠$H−κ$ and NA) are very similar, with few exceptional cases. Slight variations of those values could be a result of selection of a priori value of average $Vp$ used in the $H−κ$ method. Since, no a priori information is used in NA, the results estimated by NA are used for the results and discussion. Use of both $H−κ$ and NA makes the estimation more robust in this study. Detailed discussion is focused mainly on those cratons having sufficient stations (Pilbara, Yilgarn, Gawler, Tanzania, Kaapvaal, Zimbabwe, and Antananarivo cratons). But previous studies in these cratons have been included for comparison in crustal properties. Baranov and Morelli [45] suggested that Napier crustal thickness to be uniform of around 40 km. Here, we obtain a higher limit of ~43 km using NA (Figure 7) and ~47 km using the $H−κ$ method (Figure 9). The Moho depth of the station in Gawler craton (~46 km) is slightly higher than the highest Moho depth (~44 km) indicated by Salmon et al. [46]. They have also been considered in this study for comparative study in crustal properties. We will at first discuss the crustal thickness and $Vp/Vs$ ratio results of the six cratons (Pilbara, Yilgarn, Gawler, Tanzania, Kaapvaal, Zimbabwe, and Antananarivo cratons) obtained in this study and comparing our results with previous models of crustal studies. This is followed by discussion of crustal nature and rocks in the areas using $Vp/Vs$ and $Vs$ values. Crustal parameters have also been used to infer the degree of crustal reworking and relative density. Our results also support the link between seismic structure, evolution, and Gondwana connection.

Results estimated in our study are in agreement with the previous studies done on those areas (for example in Australia by Kennett et al. [47] and Salmon et al. [46], in Africa by Youssof et al. [48], and in Madagascar by Andriampenomanana et al. [49]). Results obtained from our study show significant variation in crustal thicknesses (~28-43 km) within and across the different Archean cratons ~28-32 km in Pilbara, ~36-42 km in Yilgarn, 36-41 km in Tanzania, 32-42 km in Kaapvaal, 35-40 km in Zimbabwe, and 36-43 km in Antananarivo (Table 1; Figure 8; Figure 10). The crustal thickness in single stations in EDC, Gawler, and Napier are ~33 km, ~46 km, and~43 km, respectively (Figure 7). The Eastern Dharwar Craton was studied extensively by Borah et al. [50] and Sarkar et al. [51]. In Napier and Gawler, along with our results, we refer to the results of Baranov and Morelli [45] and Salmon et al. [46]. These variations of crustal thicknesses are very much similar to values calculated by Durrheim and Mooney [6] in different Archean cratons (27-40 km), with a higher limit obtained in our study. The average crustal $Vs$ varies between ~3.4 and 3.8 km/s in the stations. The average crust $Vp/Vs$ ratio of all the Archean cratons varies between 1.65 and 1.78, with the exception of station BATG that has a $Vp/Vs$ value of 1.84, which indicates that the crust beneath the Archean cratons is mostly of felsic to intermediate in nature. Crustal delamination could possibly have resulted in the stabilization of the Archean cratons, leading to dominant felsic to intermediate composition of the Archean cratons [52].

The crustal thickness beneath the stations varies from 28 to 32 km with an average towards the lower limit. This suggests that the crust beneath Pilbara is thin. Our results show good correlation with Moho depth models (28-33 km) estimated in this region [46, 47, 53]. Our observations show that the $Vp/Vs$ ratio varies from 1.66 to 1.74, with an average of ~1.71. This estimation of crustal properties beneath the Pilbara in this study is consistent with the average crustal thickness of 30 km and a $Vp/Vs$ ratio of 1.70 in earlier models [54]. The average crustal shear wave velocity varies from 3.71 to 3.85 km/s. Gradational Moho is also observed beneath a few stations (Figure 7 and Supplementary) Mafic underplating was reported in previous studies ([55]; Drummond and Collins [56], which could be a possible reason for the existence of transitional Moho beneath few stations.

Previous studies have estimated the Yilgarn crust to be relatively thicker (>35 km) than Pilbara (Yuan, 2015). The crustal thickness observed in this study beneath the stations in this region is much higher than the Pilbara and varies from ~36 to 42 km. The $Vp/Vs$ ratio varies from ~1.70 to 1.76, with an average of ~1.74. This is consistent with the previous $Vp/Vs$ average model (~1.73 beneath Yilgarn) revealed by Yuan [57]. The average crustal $Vs$ varies from 3.76 to 3.8 km/s.

The crustal thickness beneath the stations in our study varies from 36 to 41 km in Tanzania. This is comparable with the average Moho depth of 39 km in the Archean terranes of Africa [58] using the gravity field model. Our results also show a good correlation with the crustal thickness values (34-40 km) estimated beneath the Tanzania craton in the Kenya Rift International seismic project [59, 60]. The crustal average $Vp/Vs$ ratio in our study varies from 1.73 to 1.77. This is consistent with the $Vp/Vs$ values (1.75-1.77) reported by Hodgson et al. [61] with a slightly lower limit. The crustal shear wave velocity beneath the stations in our study varies from 3.56 to 3.78 km/s. Gradational Moho is also observed beneath few stations (Figure 7 and Supplementary 3), indicating the presence of mafic underplating. Magmatic events such as intrusion of mafic dyke swarms were revealed in studies in the Tanzania craton [62], thus confirming the reason for the gradational Moho obtained in few stations.

The Kalahari cratons show approximately similar variations in crustal thickness (32-42 km in Kaapvaal and 35-40 km in Zimbabwe). The $Vp/Vs$ ratio varies from 1.73 to 1.75 in Kaapvaal and 1.72 to 1.77 in Zimbabwe. These values are in good correlation with crustal studies done in Africa [48]. Moho is observed clearly beneath the stations in good correlation with studies by Youssof et al. [48], except for station SA55.

The crustal thickness beneath the stations in our study varies from 36 to 42 km, with a $Vp/Vs$ ratio variation of 1.65-1.84, with an average of 1.75. This is consistent with the average $Vp/Vs$ ratio of 1.75 in Antananarivo suggested by Rindraharisaona et al. [63]. The presence of mafic lower crust in this craton [63] leads to higher value of $Vp/Vs$ (>1.8) beneath some stations. We also observe the Moho transition to be gradational in a few stations (Figure 7 and Supplementary 3). The gradational nature of Moho obtained in this study in some stations confirms the mafic underplating beneath this craton. Similar range of Moho depth (36-43 km) was also reported by Andriampenomanana et al. [49]. The crustal shear wave velocity beneath the stations in this region varies from 3.66 to 3.84 km/s. An average crustal velocity of 3.7 km/s was suggested by Andriampenomanana et al. [49]. This is consistent with our obtained value of average crustal velocity of 3.74 km/s. Pacquette et al. [33] suggested that the Neoarchean rocks in this craton were subjected to magmatism, deformation, and metamorphism that resulted in complex structured higher metamorphic grade rocks. Mafic intrusive rocks subjected to extreme conditions of metamorphosis had resulted in the formation of mafic gneisses in Madagascar cratons (Boger et al., [64]; [65, 66]), thus having a high $Vp/Vs$ ratio. This geological study is coherent with the geophysical $Vp/Vs$ ratios obtained in this study beneath the stations.

Laboratory experiments conducted by Christensen [67] indicated that rock silica percentage affects the $Vp/Vs$ ratios. The $Vp/Vs$ ratio for felsic crust having $SiO2%>66%$ is less than 1.75 (e.g., $Vp/Vs$ of felsic rock granite ~1.71); the ratio is 1.75-1.80 for intermediate crust having SiO₂%~52-66% (e.g., $Vp/Vs$ of intermediate rock diorite ~1.78), and more than 1.80 for mafic crust having SiO₂%~45-52% (e.g., $Vp/Vs$ of mafic rock gabbro is ~1.87). The presence of fluids can also affect the $Vp/Vs$ ratios, increasing it even further. This typically occurs in subducted oceanic crust [68, 69]. The $Vp/Vs$ ratio results estimated in the different cratons in this study range from 1.65 to 1.78 indicating that the nature of the crustal composition on average is felsic to intermediate. In order to check the possible rock types in the Precambrian cratons, we have averaged the laboratory data of the $Vp/Vs$ ratio and $Vs$ for the crustal layers (upper crust at pressure of 200 MPa, middle crust at pressure of 600 MPa, and lower crust at pressure of 800 MPa), with reference to Christensen [67], illustrated in Figure 11. We have then superimposed on Figure 11 the average crustal $Vp/Vs$ ratio and $Vs$ values for the stations denoted in colored circles that we have estimated using the NA modeling. The rocks inferred from the $Vp/Vs$ ratios and $Vs$ values in our study were found to be similar to the geological studies of rocks in those areas. Granitic complexes are present in the southeastern part of the Kurrana terrane, and granite greenstone units are present in the western and eastern part of the Pilbara region. The lithology of the Kalahari cratons consists of granites, greenstone belts, and gneisses, consistent with the $Vp/Vs$ ratios (~1.72). The high obtained $Vs$ values in Madagascar craton (as high as 3.8 km/s) are similar to the $Vs$ of mafic granulites (~3.8 km/s). Large belts of mafic schist and gneiss constitute the Tsaratanana complex in Antananarivo, which might result in high $Vs$ values. The high-grade structurally complex metamorphic rocks of the Neoarchean age were subjected to magmatism, deformation, and metamorphism during the Proterozoic [33]. Beneath Kaapvaal craton, biotite (tonalite) gneiss has been observed [70], leading to the observed $Vp/Vs$ and $Vs$ values as in Figure 11. Braun and Kriegsman [71] suggested the presence of felsic paragranulites in the Dharwar. Our results indicate that in Dharwar, the $Vp/Vs$ and $Vs$ values are close to paragranulite, felsic granulite, and diorite (Figure 11).

Although some undeformed ancient crust can preserve the original crustal properties, denudations alter the surface topography and the geological history of the craton. Tectonics can change the crustal properties of the ancient cratonic crust during the long evolutionary processes, such as magmatism, underplating, or delamination due to rifting or compression. Chevrot and van der Hilst [72] suggested that negative regression line between $Vp/Vs$ and Moho has the possibility to provide evidence for past tectonic activities and crustal growth such as orogenic processes or erosion. If a crustal block undergoes collision, orogeny pushes the upper crust upwards and thickens it, hence increasing the felsic component. This leads to lowering of the $Vp/Vs$ ratio. The opposite occurs in erosional process, i.e., the upper felsic part is removed, resulting in thinner mafic crust (high $Vp/Vs$ ratio). So, two opposite types of tectonic processes can result in negative correlation of $Vp/Vs$ ratio with Moho. After comparing $Vp/Vs$ with Moho (Figure 8), it was found that no concrete relation in trend of regression was observed between the $Vp/Vs$ and Moho in our study. We think that crustal reworking has a great contribution in differences in expected and observed $Vp/Vs$ values when compared to the Moho. Also, to compile the relation between $H$ and $Vp/Vs$⁠, we have also prepared Figure 12 that shows the crustal thickness (⁠$H$⁠) versus $Vp/Vs$ from previous literature for the cratons in this study superimposed with our estimations to check the correlation. We see that from having many data points of $H$ and $Vp/Vs$ from our study region, there exists no correlation between the crustal thickness and $Vp/Vs$ ratio. Relation between crustal thickness and $Vp/Vs$ is not observed globally and depends on the tectonics of the region. Cratonic regions are highly stable regions. The relations may hold for active orogenic regions [72].

Additionally, the relative crustal density can be inferred by plotting elevation versus crustal thickness, and the Airy isostatic equilibrium line assuming 3300 kg/m³ mantle density and 2760 kg/m³ average crustal density as illustrated in Figure 13. The points lying above the red line are relatively denser than isostatic equilibrium. The green triangles show the dynamic topography of the Tanzania craton. It is possibly due to the presence of a superplume beneath the craton that has led to the dynamic topography of the craton [27]. The relative crustal density is also high in most of the Tanzanian stations. For Australian cratons, the relative density of Pilbara craton is high (Figure 13), compared to Yilgarn. This also supports the vertical tectonics that operated in the region [16].

The Hadean and Archean periods were important periods during Earth’s history, and the crustal properties of the cratons belonging to the rocks of the Precambrian are important as they can give valuable clues in understanding continental structure and evolution. Current global crustal studies suggest the average crustal thickness of continents in the range of ~35-40 km. As the Precambrian cratons have survived the cycles of continental merging and rifting, it is worth investigating to see if the continental crusts of the cratonic regions follow similar range and also whether there are Moho variations with respect to the age of the Precambrian cratons as illustrated in Figure 10. We have also compiled previous estimates and included them in the plot. The rocks of the Yilgarn date back to the oldest period, Hadean. The Moho estimates are almost similar to the present-day continental estimates. The Archean eon can be divided into 4 geologic eras—the Eoarchean (4-3.6 Ga), the Paleoarchean (3.6-3.2 Ga), the Mesoarchean (3.2-2.8 Ga), and the Neoarchean (2.8-2.4 Ga). The Napier and Pilbara cratons belong to the Eoarchean, and the Moho depths fall outside the global continental average estimates. It is interesting to note the difference in Moho depth from two of the oldest Archean cratons: very shallow in Pilbara stations (~28 km-31 km) and very deep in the station of Napier (~43 km), even though both belong to the same age of 3.8 Ga. The difference in the crustal formation process in the Pilbara and the Antarctic craton could possibly have an effect on the crustal nature. It has been suggested that due to vertical tectonics operated in the Pilbara and due to partial convective overturn, the dome and keel structures were formed [16]. The Napier complex is a unique ultrahigh-temperature metamorphic complex that has ever been recorded with the highest grade of metamorphism. Various tectonothermal events, high-grade metamorphism, and plate tectonics led to Napier orogeny and formed the thick cratonic nucleus. The Paleoarchean includes the Kaapvaal and Zimbabwe cratons. Tanzania and Gawler belong to the Mesoarchean, and EDC and Antananarivo belong to the Neoarchean. Careful observation of Figure 10 shows that overall (excluding Pilbara and Gawler), there is not much variation in range of average crustal thickness in this study (~33-45 km) with respect to age of craton during the Precambrian. This range (~33-45 km) is slightly wider than the current global average ~35-40 km Moho depth. Previous studies indicate that the conditions of crustal growth have varied during the Archean and plate tectonics were dominant in the middle to late Archean. Plate tectonics control the crustal thickness to a large extent [73]. In the early Earth, the temperature was very high, and the plates were less dense and buoyant, so plate tectonics was probably not that effective [74]. Heat flow in the earth controls the rate of cooling of rocks. It takes some time for the lithosphere to cool down and the plates to become denser [75]. Neoarchean rocks were subjected to magmatism, deformation, and metamorphism during the Proterozoic [33], and underplating in the lower crust may lead to accumulation of more crustal material. But crustal melting processes like delamination can lead to melting and destruction of the lower crust. Hence, an increase in crustal thickness could not be preserved during the later part of Archean. Durrheim and Mooney [6] suggested that the mafic layer in the lower crust started appearing in the Proterozoic, with thicker crust appearing in the Proterozoic. The observations of few gradational Moho in our study beneath few stations indicate that the Archean crusts were subjected to the possibility of mafic underplating at some point, or thick crusts were present even during the Archean. In fact, Baranov and Morelli [45] reported that the central part of Antarctic Grunehogna craton shows an exceptionally thicker crustal root, which could be added later or formed initially.

Other than delamination, the role of episodic cycles of crustal evolution is also observed in the pattern of crustal thickness across each division of the Precambrian to some extent. Figure 10 shows that the range in the variation of observed present crustal thicknesses follows an interesting relatively high and low pattern alternatively with respect to age. The Hadean crust shows ~7 km variation in crustal thickness, followed by ~16-17 km in the Eoarchean and Paleoarchean. The Mesoarchean shows ~11 km Moho depth variation, increasing to ~15 km in the Neoarchean. Such a situation could take place due to alternative high and low crustal regeneration process, repeated episodic recycling, and reworking of early crust. Repeated episodes of basaltic underplating as indicated by Nelson [7] could also account for such differences in evolution of crust through time. Geodynamic models indicate that episodic cycles of large high and low melt production and crustal growth were prevalent in the early Earth [3]. Not only of the formation process, there is also an equal possibility of the other situation of later processes affecting the crustal thicknesses that we see today. The early crust after formation can go through many later processes that impact its original thickness. But the strong cratonic roots and strong lithosphere have allowed them to survive the later processes.

The cratons have survived the cycles of merging and rifting of continents. Plate tectonics have controlled the conditions of crustal formation, growth, and evolution. During Gondwana amalgamation, the collision of Archean and Proterozoic continental fragments led to Pan-African orogeny that comprises the bedrock geology of the African cratons, and this is similar to the $Vp/Vs$ and $Vs$ of the rocks in those areas. The Mozambique belt was formed as a suture between plates during the orogeny. Stern [76] suggested that it resulted from subduction and collision of east and west Gondwana. Different types of collision and also collision rate could have been an important parameter in crustal thickness variation of the cratonic areas. Crustal modification alters the composition. Geological studies indicate the presence of dominantly felsic to intermediate rocks in the cratonic areas, which is consistent with the $Vp/Vs$ ratios obtained in our study. The breakup of the supercontinent Gondwanaland led to the separation of the continents and their gradual movement in different directions. This may lead to crustal evolution and change in crustal properties. Difference in spreading and migration rate of the continents also leads to crustal delamination, and destruction of crustal roots. The present-day crustal thickness and composition were affected by plate tectonics to a large extent. Some parts of East Gondwanaland, Australia, and India moved at a faster rate than other continents. Rifting histories of Madagascar, Antarctica, and Africa also played an important role in their crustal evolution.

The Antananarivo craton in Madagascar and the African cratons—Zimbabwe and Tanzania—show a similar range of crustal thickness. Sarkar et al. [51] had studied the Dharwar craton in detail and revealed the average Moho depths of ~37 km in Dharwar craton (~34 km in Eastern Dharwar and~41 km in Western Dharwar). Sri Lanka also had drifted along with India, and an average Moho depth of ~37 km was recorded in the study by Mukherjee et al. [2]. Although the current position of Madagascar is close to Africa, Madagascar and India initially were part of a single tectonic plate during Gondwana configuration and had started rifting post Gondwana breakup. Then India separated from Madagascar and started drifting at a high speed towards Asia in the current position. The similar range of crustal thicknesses in Antananarivo, Tanzania, Zimbabwe, Dharwar, and Sri Lanka supports their Gondwana connection. The rifting history of Africa also possibly led to its crustal modification. At first, it separated along with South America from East Gondwanaland, and then later, it separated from South America giving birth to the Atlantic Ocean. The African cratons also, especially Kaapvaal, show a huge variation of thickness. The eastern part of Africa served as the connecting region between West and East Gondwana, and the cratons in this region went through crustal modifications to a large extent, resulting in varied Moho depths. Antarctica was also closer to India initially, and Southern Australia was attached to east Antarctica. Then Antarctica moved Southwards, and Australia moved eastwards at a high rate. Geochemical evidence indicates that Gawler craton was part of the Antarctic Shield prior to breakup of Gondwana [77], and the closer Moho depth in this study (~46 km in Gawler and~43 km in Napier) also supports this conclusion. Zegers et al. [28] had observed large shear zones in the Pilbara region. The Pilbara was exposed to vertical convective overturns, and the extensive extensional environment also led to crustal modifications [78]. All of this suggests a reduction in its crustal thickness, in agreement to our estimations. The crusts of the Precambrian cratons evolved to a great extent after the separation from Gondwana configuration. The present crustal properties seemed to have been affected by plate tectonics, spreading velocities and large distances covered by the continents. Australia and India covered a large distance after splitting up completely from Gondwana.

This work analyses crustal properties across the ancient Precambrian cratonic regions that previously made up the Gondwana supercontinent. Using the teleseismic analysis method of the receiver function analysis, we have analysed and compiled the crustal thickness, average $Vp/Vs$⁠, ratios and crustal $Vs$⁠. These are then compared across cratonic locations, and we attempt to infer details of ancient crustal formation processes and subsequent crustal remodeling.

• (i)

  On comparing the crustal thickness with the age of the cratons, it was found that average bounds of crustal thickness were similar between Hadean, early Archean, and late Archean cratons (~33-45 km). The cratons have a dynamic history of migration and separation after Gondwana breakup with crustal modifications. Even if plate tectonics was dominant in the middle to late Archean, difference in spreading and migration rate of the continents led to crustal delamination, and destruction of crustal roots. Comparison of the crustal properties was also found to support their Gondwana connection

• (ii)

  There are intracratonic variations of Archean crustal thicknesses beneath the studied stations of eastern Gondwana: ~28-32 km in Pilbara; ~36-42 km in Yilgarn; ~46 km in Gawler; 36-41 km in Tanzania; 32-42 km in Kaapvaal; ~35-40 km in Zimbabwe; and ~36-43 km in Antananarivo. It also shows nonuniform crustal thickness variation of the Archean crusts. The stations in this study show the thickness ranges from less than 30 km in parts of Pilbara to more than 40 km beneath parts of Yilgarn, Tanzania, Kaapvaal, Zimbabwe, and Antananarivo. Thick crusts of early Archean age suggest the possibility of high rate of crustal formation or crustal modification, due to which materials could be added later leading to crustal thickening. $Vp/Vs$ ratios mostly in the range ~1.65-1.78 indicate felsic to intermediate nature of the Archean cratons, with evidence of gradational Moho beneath few stations. The rocks inferred from the $Vp/Vs$ ratios and $Vs$ values in our study and silica percentages were found to be similar to the global laboratory experiments and geological studies of rocks in those areas

• (iii)

  Other than delamination, the role of episodic cycles of crustal growth and evolution is also observed in the pattern of crustal thickness observed at present stage in each division of the Precambrian. The Hadean crust shows ~7 km variation in crustal thickness, followed by ~16-17 km in the Eoarchean and Paleoarchean. The Mesoarchean shows ~11 km Moho depth variation, increasing to ~15 km in the Neoarchean. A number of factors could account for such oscillating differences: operation of plate tectonics through time, alternative high and low crustal regeneration process, repeated episodic recycling, and reworking of crust. Investigation of a new dense network of stations in the Precambrian cratons and integrating numerous geological and geophysical detailed observations in future would help in understanding early Earth in more detail.

Seismic waveform data included in this study can be obtained from the following: Incorporated Research Institutions for Seismology (IRIS) Data Management Centre (DMC), Institut de Physique du Globe de Paris (IPGP), and Réseau sismologique et géodésique français (RESIF) data centers. The networks are as follows: (i) G (doi:10.18715/GEOSCOPE.G), (ii) AU (https://www.fdsn.org/networks/detail/AU/), (iii) XD (doi:10.7914/SN/XD_1994), (iv) XA (doi:10.7914/SN/XA_1997), (v) II (doi:10.7914/SN/II), and (vi) XV (doi:10.7914/SN/XV_2011). Data citations for the above networks are given as follows: (i) Network G (network name: GEOSCOPE). Institut De Physique Du Globe De Paris (IPGP), & Ecole Et Observatoire Des Sciences De La Terre De Strasbourg (EOST) (1982). GEOSCOPE, French Global Network of broadband seismic stations. Institut de Physique du Globe de Paris (IPGP), Université de Paris (doi:10.18715/GEOSCOPE.G). (ii) Network AU (network name: Australian National Seismograph Network ANSN). Operated by Geoscience Australia (https://www.fdsn.org/networks/detail/AU/). (iii) Network XD (network name: Seismic investigations of the Lithospheric Structure of the Tanzanian Craton (Tanzania)). Tom Owens and Andy Nyblade (1994). Seismic investigations of the Lithospheric Structure of the Tanzanian Craton (data set). International Federation of Digital Seismograph Networks (doi:10.7914/SN/XD_1994). (iv) Network XA (network name: Anatomy of an Archean Craton, South Africa, a Multidisciplinary Experiment across the Kaapvaal Craton (SASEK)). Paul Silver (1997). Anatomy of an Archean Craton, South Africa, a Multidisciplinary Experiment across the Kaapvaal Craton (data set). International Federation of Digital Seismograph Networks (doi:10.7914/SN/XA_1997). (v) Network II (network name: Global Seismograph Network- IRIS/IDA (GSN)). Scripps Institution of Oceanography (1986). Global Seismograph Network-IRIS/IDA (data set). International Federation of Digital Seismograph Networks (doi:10.7914/SN/II). (vi) Network XV (network name: INVESTIGATION OF SOURCES OF INTRAPLATE VOLCANISM USING PASSCAL BROADBAND INSTRUMENTS IN MADAGASCAR, THE COMORES, AND MOZAMBIQUE (MACOMO)). Michael Wysession, Douglas Wiens, & Andy Nyblade (2011). INVESTIGATION OF SOURCES OF INTRAPLATE VOLCANISM USING PASSCAL BROADBAND INSTRUMENTS IN MADAGASCAR, THE COMORES, AND MOZAMBIQUE (data set). International Federation of Digital Seismograph Networks (doi:10.7914/SN/XV_2011).

The authors declare that they have no conflicts of interest.

Seismic waveform data has been obtained from the following: Incorporated Research Institutions for Seismology (IRIS) Data Management Centre (DMC), Institut de Physique du Globe de Paris (IPGP), and Réseau sismologique et géodésique français (RESIF) data centers. The networks are as follows: (i) G (doi:10.18715/GEOSCOPE.G), (ii) AU (https://www.fdsn.org/networks/detail/AU/), (iii) XD (doi:10.7914/SN/XD_1994), (iv) XA (doi:10.7914/SN/XA_1997), (v) II (doi:10.7914/SN/II), and (vi) XV (doi:10.7914/SN/XV_2011). Data citation is included after this section. Most of the data processing was done using Seismic Analysis Code (SAC), and figures were made using Generic Mapping Tools (GMT). Author P. Mukherjee is grateful to the INSPIRE fellowship, Department of Science and Technology (DST), India, and MEXT fellowship, Japan Science and Technology (JST), for supporting the study. K. Borah gratefully acknowledges the financial support by the Science and Engineering Research Board (SERB), India (Project No. MTR/2019/001260 and CRG/2020/003726), and Academic Research Funding (ARF) from IISER Kolkata. P. Mukherjee thanks others from IISER Kolkata Earth Science Department who have helped in this study. P. Mukherjee thanks Dr. Y. Ito for guidance and research budget.