The computational costs of sensitivity analysis for a considerably large simulation is attributed to two main phenomena. Firstly, as it's mentioned in the body of the question, due to a large number of parameters, even a small selection of samples from each parameter would lead to an astronomical number of runs. Secondly, a small sampling number (two-level experiment) is never enough to even roughly understand the behavior of a system with regards to its parameters.
After spending a couple of days, this is what I understood. There are usually two steps taken to address the raised concerns. First, instead of using a full factorial design, one can use the fractional factorial design to reduce the number of runs to address the first concern. This approach assumes that most of the information provided by the full factorial design is redundant and can be ignored accordingly. Read more on Wikipedia. Using the fractional factorial design with a two-level factor, one can conduct a primary study on parameters and drop those that are insignificant. This would result in a much smaller parameter set that can be more extensively investigated using other methods (mentioned below). For the sake of completeness, there are also other methods in this scope to reduce the complexity of full factorial design such as Plackett-Burman designs, Cotter designs, and mixed-level designs.
In the second step, in order to further understand the system, one can use methods such as those classified as Space-filling design (e.g. Latin hypercube sampling-LHD) to investigate each parameter's influence on the system in more detail. Different methods use different sampling approaches based on the nature of a problem. I find this open article ("Sensitivity analysis by design of experiments") by Schepdael et al 2016 useful in this regard.
One more point. If running your model takes a substantial amount of time and executing thousands of runs based on the aforementioned methods are not feasible, check out Approximate Approximate Bayesian Calculation. This method develops a statistical model based on the limited result sets from your real model an executes more extensive runs on the statistical model instead. I think this is the main article elaborating this method: "AABC: Approximate approximate Bayesian computation for inference in population-genetic models" by Buzbas et al 2015.
Please correct me or let me know your opinions in the comments.