In the framework of risk assessment in nuclear accident analysis, best-estimate computer codes associated with probabilistic modeling of uncertain input variables are used to estimate safety margins. Often, a first step in such uncertainty quantification studies is to identify the critical configurations (or penalizing, in the sense of a prescribed safety margin) of several input parameters (called scenario inputs) under the uncertainty of the other input parameters. However, the large CPU-time cost of most of the computer codes used in nuclear engineering, as the ones related to thermal-hydraulic accident scenario simulations, involves developing highly efficient strategies. This work focuses on machine learning algorithms by way of a metamodel-based approach (i.e., a mathematical model that is fitted on a small sample of simulations). To achieve it with a very large number of inputs, a specific and original methodology called Identification of penalizing Configurations using SCREening And Metamodel (ICSCREAM) is proposed. The screening of influential inputs is based on an advanced global sensitivity analysis tool (Hilbert-Schmidt Independence Criterion importance measures). A Gaussian process metamodel is then sequentially built and used to estimate within a Bayesian framework the conditional probabilities of exceeding a high-level threshold according to the scenario inputs. The efficiency of this methodology is illustrated with two high-dimensional (around a hundred inputs) thermal-hydraulic industrial cases simulating an accident of primary coolant loss in a pressurized water reactor. For both use cases, the study focuses on the peak cladding temperature (PCT), and critical configurations are defined by exceeding the 90%-quantile of the PCT. In both cases, using only around one thousand code simulations, the ICSCREAM methodology allows one to estimate the impact of the scenario inputs and their critical areas of values.