Bayesian inversion is commonly applied to quantify uncertainty of hydrological variables. However, the focus in Bayesian inversion is more on spatial hydrological properties instead of hyperparameters or global/non-gridded variables. In this paper, we present a hierarchical Bayesian framework to quantify uncertainty of both global and spatial variables. We estimate first the posterior of global variables and then hierarchically estimate the posterior of the spatial field. We propose a machine learning-based inversion method to estimate the joint distribution of data and global variables directly without introducing a statistical likelihood. We also propose a new local dimension reduction method: local principal component analysis (local PCA) to update large-scale spatial fields with local data more efficiently. We illustrate the hierarchical Bayesian formulation with two case studies: one with a linear forward model (volume averaging inversion) and a second with a non-linear forward model (pumping tests). Results show that quantifying global variables uncertainty is critical for assessing uncertainty on predictions. We show how the local PCA approach accelerates the inversion process. Furthermore, we provide an open-source Python package on the hierarchical Bayesian framework including two case studies.