Abstract: This study aims to systematically evaluate the performance of multiple medium-to-long-term runoff forecasting models at nine major river control stations in Sichuan Province, China, and explore ensemble methods to enhance forecast accuracy and robustness. Six models were selected: a linear model (Generalized Linear Model, GLM), machine learning models (Random Forest, RF; Boosted Regression Trees, BRT; Cubist), and neural network models (Long Short-Term Memory, LSTM; Back Propagation, BP). Utilizing observed runoff and precipitation data (1980-2024), along with 88 atmospheric circulation indices and 26 sea surface temperature indices, model parameters were calibrated via Bayesian optimization combined with 10-fold cross-validation (calibration period: 1980-2008; validation period: 2009-2024) to assess generalization capability. High-performing models were integrated using Bayesian Model Averaging (BMA), Quantile Model Averaging, and Ensemble Selection. Feature importance analysis was conducted to identify key drivers of runoff. Results demonstrated that LSTM, RF, BRT, and Cubist models significantly outperformed the traditional GLM and shallow BP model across multiple metrics. During calibration, RF achieved optimal performance (median Nash-Sutcliffe Efficiency, NSE = 0.91; median Mean Absolute Percentage Error, MAPE = 13.7%). In validation, LSTM (median NSE decreased from 0.79 during calibration to 0.63) and BRT (median MAPE increased from 21.7% to 25.2%) exhibited relatively stronger stability. Feature importance analysis revealed the Tibetan Plateau Index-1 and the Asian Polar Vortex Intensity Index as the dominant controlling factors for runoff at the Wutongqiao Station on the Min River. The impact of the same climatic index with different time lags on runoff could be contrasting, linked to its nonlinear interaction with seasonal adjustments in moisture transport pathways. The BMA ensemble model showed relative performance, matching the best single model (RF) in validation period median NSE (0.67) and ranking second only to BRT in median MAPE (25.2%), effectively constraining hydrological uncertainty. This study provides optimization solutions for runoff simulation in complex basins through multiple model comparison and integration, and the hybrid modeling strategy enhances the robustness of predictions, demonstrating significant practical significance.