Passive Safety Systems (PSSs) are increasingly employed in advanced Nuclear Power Plants (NPPs) since they are considered, in general, more reliable than active systems. Their safety performance is evaluated through computationally expensive Thermal-Hydraulic (T-H) simulations models and the identification of the operational conditions which lead to unsafe conditions (the so-called Critical failure Regions, CRs) may be challenging. In this view, the computationally expensive T-H models simulating the PSSs behavior can be replaced by fast-running surrogate models (also called metamodels), coupled with adaptive sampling techniques for speeding up the exploration and efficiently focusing the analysis on the most interesting regions of the domain, i.e., the CRs boundary (limit surface).