## 3 Numerical Framework
### Shock Capturing for Material Interfaces  
The VDI flow model allows dealing with incompressible flow problems involving immiscible fluids of different densities, e.g., water and air. Different fluid components are discriminated through the density, which acts as a phase variable, and the interface between components is represented as a sharp variation of the density field. It is well known that high-order discretizations are affected by Gibbs-type phenomena in the presence of excessively sharp gradients. This might endow a lack of monotonicity at interfaces, which must be properly controlled to avoid stability issues. In this work, we implemented a strategy based on the introduction of artificial viscosity inspired by the approach proposed by [PERSON] and [PERSON] (1992) in the context of compressible flows featuring shocks. In the context of the VDI flow model, spurious density oscillations might occur at contact discontinuities. Accordingly, the artificial diffusive term is introduced solely within the continuity equation  
\[\frac{\partial\rho}{\partial t}+\frac{\partial}{\partial x_{j}}(\rho u_{j})= \frac{\partial}{\partial x_{j}}\left(v_{e}\frac{\partial\rho}{\partial x_{j} }\right).\]  
Nevertheless, in order to ensure that the momentum equation of (1) is actually solved, we propose the following conservative form of the VDI flow model with artificial viscosity:\[\frac{\partial u_{j}}{\partial x_{j}} =0,\] \[\frac{\partial}{\partial t}(\rho u_{i})+\frac{\partial}{\partial x_{j} }(\rho u_{i}u_{j}) =-\frac{\partial p}{\partial x_{i}}+\frac{\partial\tau_{ij}}{ \partial x_{j}}+\rho g_{i}+u_{i}\frac{\partial}{\partial x_{j}}\left(v_{e} \frac{\partial\rho}{\partial x_{j}}\right), \tag{14}\] \[\frac{\partial\rho}{\partial t}+\frac{\partial}{\partial x_{j}}( \rho u_{j}) =\frac{\partial}{\partial x_{j}}\left(v_{e}\frac{\partial\rho}{ \partial x_{j}}\right),\]  
with \(i,j=1,\ldots,d\) and \(v_{e}\) is the artificial mass diffusion parameter. In order to ensure that artificial mass diffusion is introduced only at interfaces, we rely on a smoothness sensor based on the modal decay of the working density \(\widetilde{\rho}\) expansion. The artificial diffusion \(v_{e}\) is, thus, defined according to [34] as  
\[v_{e}= \begin{cases}\
u_{e0}&\text{if}\quad s<s_{0}-\Delta s\\ \frac{v_{e0}}{2}\Big{[}1+\sin\left(\frac{\pi}{2}\frac{s_{0}-s}{\Delta s} \right)\Big{]}&\text{if}\quad s_{0}-\Delta s\leq s\leq s_{0}+\Delta s\\ 0&\text{otherwise}.\end{cases} \]  
where the smoothness sensor \(s\) is defined according to the _skyline pessimization_ procedure proposed by [PERSON] et al. [35]. According to [35], \(s\gg 1\) when the field is smooth and \(s\approx 1\) when the density field exhibits abrupt variations. The user-defined parameters \(s_{0}\) and \(\Delta s\) define the threshold triggering the introduction of artificial viscosity as well as the actual amount of artificial viscosity introduced, which varies from \(0\) to \(v_{e}\). In this work, we set \(s_{0}=1\), \(\Delta s=0.75\) and \(v_{e0}=h/k\), with \(h\) the mesh element size and \(k\) the DG polynomial degree. This choice aims to introduce artificial diffusion that decreases with higher polynomial degrees and smaller elements, allowing to possibly capture discontinuities as thin layers within a single mesh element. A sufficiently high polynomial degree, namely \(k\geq 4\), is required to ensure that the modal decay provides a reliable smoothness sensor. The two additional terms related to artificial diffusion appearing in (14), one within the continuity equation and one within the momentum equation, are discretized in space following the idea proposed by [PERSON] et al. [36]. In particular, the additional term within the continuity equation is discretized by multiplying by an arbitrary test function, integrating by parts and neglecting the numerical flux contributions on the inter-element boundaries.