## 1 Introduction  
A thorough analysis of how a system's parameters influence its response is crucial for optimal and robust product design, which reduces time and costs in the late stages of product development. Running high-fidelity simulations of the physical system, however, requires a large amount of computational resources; thus, extensive exploration over the design space is prohibitive. Because of this, surrogate models have been used in different domains as efficient tools for decision making and risk management [1, 2, 3]. Surrogates or metamodels are data-driven methods that extract and learn information from input-output pairs obtained from a physics-based simulation to create a proxy with short evaluation time and little loss of accuracy. Surrogate models have been applied to the vibroacoustic domain as noise, vibration and harshness (NVH) performance arises as a key indicator of customer satisfaction and vibroacoustic simulations are computationally costly due to, for instance, the complex behavior involved in fluid-structure interactions. In many studies, the system response is approximated by a polynomial using the response surface methodology (RSM) [4, 5, 6]. Second-order polynomials are usually used in RSM for their sample efficiency and interpretability, but they are unable to capture arbitrary nonlinearities in the system response. On the other hand, higher-order polynomials are rarely employed due to the costly training [7]. Therefore, machine learning (ML)-based surrogates have been used in the vibroacoustic analysis for their capability of learning complex data relations and nonlinear functions and applied in optimization [8, 9, 10], uncertainty analyses [11, 12, 13] and knowledge discovery from measured databases [14]. A surrogate model using Gaussian process regressor (GPR) was implemented in [9] to optimize the sound transmission loss (STL) of an intake system. [PERSON] et al. [15] trained a GPR-based surrogate model to predict the increase of the STL response at resonance frequencies of porous acoustic metamaterials. The results showed that better accuracy can be achieved when a more feature-rich STL model is used. [PERSON] et al. [10] performed the optimization of the band-gap of acoustic metamaterials using a radial basis function (RBF)-based surrogate. Surrogates for crash and NVH response of a body-in-white model using RBF and neural network (NN)-based algorithms were designed in [8; 16], respectively. [PERSON] et al. [17] investigated various ML-based surrogate models to predict the noise radiated by electric motors and concluded that global models encompassing structural, electromagnetic, and acoustic physical domains performed better than local surrogates. NN-based surrogate models have also accurately replaced boundary element methods to predict interior vehicle sound-pressure level [18]. The growing trend in vibroacoustics of relying on surrogates based on high-fidelity simulations not just enables the search for optimal and reliable designs, but also paves the way to the construction of digital twins for real-time vibroacoustic applications, such as online monitoring and control [19; 20]. However, the reliability of these solutions depends on the accuracy of the surrogate models. In particular, in vibroacoustics, the typical non-smooth functions pose difficulty to the construction of accurate ML-based surrogates as ML methods assume local smoothness of data [21; 22]. Therefore, to ensure that further surrogate-based developments on vibroacoustics are well-founded, it is imperative to develop scientific-based guidelines to construct physically consistent and accurate ML-based surrogates considering the particulars of the domain. Thus, the objective of this paper is to identify ML methods capable of modeling non-smooth functions, as well as techniques to incorporate domain knowledge into the surrogate.