Within ten nearby (d <  450 pc) Gould belt molecular clouds we evaluate statistically the relative orientation between the magnetic field projected on the plane of sky, inferred from the polarized thermal emission of Galactic dust observed by Planck at 353 GHz, and the gas column density structures, quantified by the gradient of the column density, NH. The selected regions, covering several degrees in size, are analysed at an effective angular resolution of 10′ FWHM, thus sampling physical scales from 0.4 to 40 pc in the nearest cloud. The column densities in the selected regions range from NH≈ 1021 to1023 cm-2, and hence they correspond to the bulk of the molecular clouds. The relative orientation is evaluated pixel by pixel and analysed in bins of column density using the novel statistical tool called “histogram of relative orientations”. Throughout this study, we assume that the polarized emission observed by Planck at 353 GHz is representative of the projected morphology of the magnetic field in each region, i.e., we assume a constant dust grain alignment efficiency, independent of the local environment. Within most clouds we find that the relative orientation changes progressively with increasing NH, from mostly parallel or having no preferred orientation to mostly perpendicular. In simulations of magnetohydrodynamic turbulence in molecular clouds this trend in relative orientation is a signature of Alfvénic or sub-Alfvénic turbulence, implying that the magnetic field is significant for the gas dynamics at the scales probed by Planck. We compare the deduced magnetic field strength with estimates we obtain from other methods and discuss the implications of the Planck observations for the general picture of molecular cloud formation and evolution.