This study applies Bayesian inversion to receiver functions (RFs) to estimate local shear‐wave velocity (VS) structure of the Juan de Fuca (JdF) plate beneath the northern Cascadia subduction zone (CSZ) offshore and onshore Vancouver Island, British Columbia, Canada. We use passive seismic data recorded on NEPTUNE (NorthEast Pacific Time‐series Undersea Networked Experiments) Canada ocean‐bottom seismometers (OBSs), on temporary autonomous KECK Foundation OBSs, and on two land‐based seismometers on Vancouver Island that are part of the Canadian National Seismograph Network (CNSN). Three‐component, broadband recordings of large (), distant (30°–100°) earthquakes are used to compute RFs dominated by locally generated P‐to‐S converted waves. These are subsequently inverted using a nonlinear Bayesian approach that yields optimal profiles of VS, VP (compressional‐wave velocity), and strike and dip angles, as well as quantitative uncertainty estimates for these parameters. The introduction of NEPTUNE Canada helps fill a gap in offshore seismic monitoring. Results from OBS stations indicate a thin oceanic crust at the JdF Ridge which thickens to ∼10 km at the continental slope where sediment thickness also increases to ∼5 km. At OZB, a coastal station, a 6–8 km thick, two‐part low‐velocity zone (LVZ) is imaged at 19 km depth. An LVZ of similar thickness is also observed 34 km beneath PGC, a south‐central Vancouver Island station. The thickness of the LVZ imaged at these two land‐based stations indicates that the oceanic sediments are not subducted but are scraped off the JdF plate and accreted to the North American plate. Determining these VS models at various stages of the CSZ provides a more detailed image of the subducting plate, and therefore contributes valuable new information useful for seismic‐hazard analysis.