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 H-2 clusters on fullerenes, following a first-principles approach in which the interaction is quantitatively evaluated for the C-20 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 H-2 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 C-60-H-2. 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 C-60 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.