New, full-dimensional potential energy surfaces (PESs), obtained using precise least-squares fitting of high-level electronic energy databases, are reported for intrinsic H2(H2O) two-body and H2(H2O)2 three-body potentials. The database for H2(H2O) consists of approximately 44 000 energies at the coupled cluster singles and doubles plus perturbative triples (CCSD(T))-F12a/haQZ (aug-cc-pVQZ for O and cc-pVQZ for H) level of theory, while the database for the three-body interaction consists of more than 36 000 energies at the CCSD(T)-F12a/haTZ (aug-cc-pVTZ for O, cc-pVTZ for H) level of theory. Two precise potentials are based on the invariant-polynomial technique and are compared to computationally faster ones obtained via "purified" symmetrization. All fits use reduced permutational symmetry appropriate for these non-covalent interactions. These intrinsic potentials are employed together with existing ones for H2, H2O, and (H2O)2, to obtain full PESs for H2(H2O) and H2(H2O)2. Properties of these full PESs are presented, including a diffusion Monte Carlo calculation of the zero-point energy and wavefunction, and dissociation energy of the H2(H2O) dimer. These PESs together with an existing one for water clusters are used in a many-body representation of the PES of hydrogen clathrate hydrates, illustrated for H2@(H2O)20. An analysis of this hydrate is presented, including the electronic dissociation energy to remove H2 from the calculated equilibrium structure.