Abstract The cosmic microwave background (CMB) radiation B-mode polarization signal contains the unique signature of primordial metric perturbations produced during the inflation. The separation of the weak CMB B-mode signal from strong foreground contamination in observed maps is a complex task, and proposed new-generation low-noise satellite missions compete with the weak signal level of this gravitational background. In this paper, for the first time, we employ a foreground model-independent internal linear combination (ILC) method to reconstruct the CMB B-mode signal using simulated observations over large angular scales of the sky of six frequency bands of the future-generation CMB mission Probe of Inflation and Cosmic Origins (PICO). We estimate the joint CMB B-mode posterior density following the interleaving Gibbs steps of the B-mode angular power spectrum and cleaned map samples using the ILC method. We extend and improve the earlier reported Bayesian ILC method to analyze weak CMB B-mode reconstruction by introducing noise-bias corrections at two stages during the ILC weight estimation. By performing 200 Monte Carlo simulations of the Bayesian ILC method, we find that our method can reconstruct the CMB signals and the joint posterior density accurately over large angular scales of the sky. We estimate the Blackwell–Rao statistics of the marginal density of the CMB B-mode angular power spectrum and use them to estimate the joint density of the tensor-to-scalar ratio r and a lensing power spectrum amplitude A lens . Using 200 Monte Carlo simulations of the delensing approach, we find that our method can achieve an unbiased detection of the primordial gravitational-wave signal r with more than 8σ significance for levels of r ≥ 0.01 for input frequency maps with realistic lensing and foregrounds with constant spectral indices.