Numerical simulation is important for the risk assessment of subaerial landslide-induced tsunami (SLIT) hazards; however, their predictive performances are thus far limited as these works fail to consider the involved uncertainties. We propose a novel probabilistic risk assessment framework for quantifying several uncertainties, such as the model parameters and bias. We successfully couple the numerical simulation and Bayesian-based back analysis method. The wave height and run-up are captured by a three-dimensional granular-flow landslide and renormalization group turbulence model in FLOW3D; furthermore, the relevant parameter information is calibrated using a Markov chain Monte Carlo algorithm with the differential evolution adaptive metropolis. In this framework, which incorporates the prior information of model parameters, and field observations, as well as numerical model bias, calibrated parameters exhibit a probability distribution with a small standard deviation, facilitating improved risk prediction capabilities; In addition, the optimized adaptive Kriging-based surrogate model with multiple outputs is established to reduce the computation cost of the framework with acceptable accuracy. The framework's outputs, including a spatial probability that tsunamis inundate elements at risk, and a hazard zonation map, contribute to the evaluation of SLIT hazards and prioritize mitigation measures. Our work introduces a paradigm for the application of the framework, demonstrated through studies of two landslides in Wu Gorge. Without loss of generality, the framework offers a novel perspective for the quantitative assessment of SLIT hazards.