The interaction between hydrogen and carbonaceous nanostructures is of fundamental interest in various areas of physical chemistry. In this contribution we have revisited the physisorption of hydrogen molecules and H2 clusters on fullerenes, following a first-principles approach in which the interaction is quantitatively evaluated for the C20 system using high-level electronic structure methods. Relative to coupled cluster data at the level of single, double, and perturbative triple excitations taken as a benchmark, the results for rotationally averaged physisorbed H2 show a good performance of MP2 variants and symmetry-adapted perturbation theory, but significant deviations and basis set convergence issues are found for dispersion-corrected density functional theory. These electronic structure data are fitted to produce effective coarse-grained potentials for use in larger systems such as C60-H2. Using path-integral molecular dynamics, the potentials are also applied to parahydrogen clusters solvated around fullerenes, across the regime where the first solvation shell becomes complete and as a function of increasing temperature. For C60 our findings indicate a sensible dependence of the critical solvation size on the underlying potential. As the temperature is increased, a competition is found between the surface and radial expansions of the solvation shell, with one molecule popping away at intermediate temperatures but getting reinserted at even higher temperatures.