Continental crust and the underlying lithospheric mantle make up the continental lithosphere of the Earth. Our understanding of its structure and composition is limited by the inaccessibility of Earth’s deep interior. Seismic imaging utilizing complementary seismic data provides unique constraints on the present-day structure of continental lithosphere. However, while recent efforts have improved the resolution of seismic images, the quantification of uncertainties remains challenging due to the non-linearity and the non-uniqueness of the geophysical inverse problem. To gain insights into the composition, formation, and evolution of the continental lithosphere, an interdisciplinary approach that incorporates seismological, geodynamical, and geochemical contributions is needed. In this dissertation, I implement a model-space search approach – transdimensional Bayesian inversion – to explore seismological constraints of continental lithosphere. I utilize seismic observables including Rayleigh and Love wave dispersion, Rayleigh wave ZH ratio, and Ps receiver functions to invert for shear velocity (Vs), compressional velocity (Vp), density (ρ), and radial anisotropy (ξ) profiles of lithospheric structure. I begin by systematically investigating the effects of parameterization choices on inversion results using synthetic tests. Then, I proceed to tackle several technical challenges regarding the accurate retrieval of multi-parameter velocity structures from large seismic arrays in the presence of sediment layers. Finally, I apply these techniques to create a shear velocity model (TBI-NGP) of the lithosphere across the Northern Great Plains of the United States using data from the EarthScope Transportable Array. This probabilistic seismic model enables statistical assessment of the elastic properties in an Archean craton and Paleoproterozoic orogen. Subsequently, I incorporate seismic constraints with geophysical and geochemical measurements to infer the composition of continental crust in the region.