The Superior craton preserves a geological history spanning >3 billion years, making it an exceptional natural laboratory to explore the evolution of the Earth’s ancient crust. Despite its geological importance, our understanding of the crustal structure in this region remains relatively sparse, mainly due to limited data sets and basic analysis methods used in past studies. We use earthquake data recorded at 16 Canadian seismograph stations over the period 2003-2023. P-wave receiver functions (PRF) are computed using an iterative time-domain deconvolution method, after applying a 1 Hz low-pass Gaussian filter. H-k stacking is employed to obtain initial estimations of the crustal thickness and the Vp/Vs ratio of the crust. Estimated crustal thickness varies from ∼33 km to ∼45 km through most of the province and bulk crustal composition, as inferred from Vp/Vs ratios, ranges from felsic to intermediate. To infer detailed crustal shear-wave velocity, we perform a joint inversion of the PRF with Rayleigh wave phase-velocity dispersion data within several back-azimuthal bins. This approach highlights structural variations for different directions around each station, providing valuable insights into the 3D crustal architecture. We conduct harmonic decomposition for stations with a high azimuthal coverage to investigate possible sources of anisotropy. Significant variations in crustal characteristics are observed across the region, from relatively simple structure and a sharp Moho to a layered crust and more gradual crust-mantle transition. We explore the variability in crustal structure and anisotropy in the context of the assembly and tectonic evolution of the Superior craton.