Inverse Uncertainty Quantification (IUQ) has gained increasing attention in the field of nuclear engineering, especially nuclear thermal-hydraulics (TH), where it serves as an important tool for quantifying the uncertainties in the physical model parameters (PMPs) while making the model predictions consistent with the experimental data. In this paper, we present an extension to an existing Bayesian inference-based IUQ methodology by employing a hierarchical Bayesian model and variational inference (VI), and apply this novel framework to a real-world nuclear TH scenario. The proposed approach leverages a hierarchical model to encapsulate group-level behaviors inherent to the PMPs, thereby mitigating existing challenges posed by the high variability of PMPs under diverse experimental conditions and the potential overfitting issues due to unknown model discrepancies or outliers. To accommodate computational scalability and efficiency, we utilize VI to enable the framework to be used in applications with a large number of variables or datasets. The efficacy of the proposed method is evaluated against a previous study where a No-U-Turn-Sampler was used in a Bayesian hierarchical model. We illustrate the performance comparisons of the proposed framework through a synthetic data example and an applied case in nuclear TH. Our findings reveal that the presented approach not only delivers accurate and efficient IUQ without the need for manual tuning, but also offers a promising way for scaling to larger, more complex nuclear TH experimental datasets.