Land surface models (LSMs) implemented in numerical weather prediction (NWP) models use several parameters to suitably describe the surface and its interaction with the atmosphere, whose determination is often affected by many uncertainties, strongly influencing simulation results. However, the sensitivity of meteorological model results to these parameters has not yet been studied systematically, especially in complex terrain, where uncertainty is expected to be even larger. This work aims at identifying critical LSM parameters influencing the results of NWP models, focusing in particular on the simulation of thermally-driven circulations over complex terrain. While previous sensitivity analyses employed offline LSM simulations to evaluate the sensitivity to surface parameters, this study adopts an online coupled approach, utilizing the Noah-MP LSM within the Weather Research and Forecasting (WRF) model. To overcome computational constraints, a novel tool, Machine Learning-based Automated Multi-method Parameter Sensitivity and Importance analysis Tool (ML-AMPSIT), is developed and tested. This tool allows users to explore the sensitivity of the results to model parameters using supervised machine learning regression algorithms, including Random Forest, CART, XGBoost, SVM, LASSO, Gaussian Process Regression, and Bayesian Ridge Regression. These algorithms serve as fast surrogate models, greatly accelerating sensitivity analyses while maintaining a high level of accuracy. The versatility and effectiveness of ML-AMPSIT enable the fast implementation of advanced sensitivity methods, such as the Sobol method, overcoming the computational limitations encountered in expensive models like WRF. The suitability of this tool to assess model’s sensitivity to the variation of specific parameters is first tested in an idealized sea breeze case study where six surface parameters are varied. Then, the analysis focuses on the evaluation of the sensitivity to surface parameters in the simulation of thermally-driven circulations in a mountain valley. Specifically, an idealized three-dimensional topography consisting of a valley-plain system is adopted, analyzing a complete diurnal cycle of valley and slope winds. The analysis focuses on all the key surface parameters governing the interactions between NoahMP and WRF. The proposed approach, novel in the context of LSM-NWP model coupling, draws from established applications of machine learning in various Earth science disciplines, underscoring its potential to improve the estimation of parameter sensitivities in NWP models.