## 5 Conclusions  
In this work, we investigated the numerical performance of an implicit high-order discontinuous Galerkin solver for variable density incompressible flow problems. Within the VDI model, the density is treated as a purely advected property and, thus, can be used to track possibly multiple (more than two) components. Moreover, the interface between fluids is captured in a diffuse fashion as a smooth, but sharp, variation in the density field, thus not requiring any geometrical reconstruction. A Riemann solver, namely the SDRS or the ERS proposed in [22], is used for the definition of the Godunov numerical fluxes. Both solvers are based on an artificial compressibility modification of the inviscid incompressible set of equations, which introduces a pressure-velocity coupling at the discrete level. As the DG-VDI formulation aims to deal with multi-component flow problems featuring high-density ratios (water-air), a numerical strategy ensuring density positivity at a discrete level  
Figure 8: Impinging droplet—density field for setup A (**left column**), B (**middle column**) and C (**right column**) at non-dimensional times 0.45, 0.50, 0.55, and 0.70 from top to bottom. Color scale denotes water with orange and air with blue (see Figure 7 for reference). was developed. The density is specified as a function of a different variable, the so-called _working density_\(\widetilde{\rho}\), which is part of the set of unknowns of the DG formulation. The change of density variable \(\rho=f(\widetilde{\rho})\) is suitably designed to enforce the desired property on the density. As the main drawback, the variable transformation introduces strong nonlinearities that must be accurately integrated both in space and time to guarantee mass conservation. In order to control spurious oscillations arising from the high-order polynomial description of contact discontinuities, i.e., material interfaces, an artificial mass diffusion term is introduced within the mass and momentum equations. The proposed implementation of such a term only modifies the higher modes of the DG discretization within each mesh element leaving untouched the first mode, namely the mean value of physical quantities. Time integration is performed with high-order single-step multi-stage ESDIRK schemes with optimal stability properties, which prevent Courant-related time step restrictions. The assessment of the approach in terms of space, time, and mass conservation accuracy was performed for the Kovasznay test case, the Guermond-Quaratpelle manufactured solution, and sine and square density waves. Finally, the solver was successfully used for the simulation of the impingement of a 2D liquid (water) droplet on a thin film immersed in air. Several grids with different initial interface thicknesses, i.e., different initial density gradient magnitudes at the interfaces, were considered and compared. Future work will further assess the proposed DG approach to three-dimensional flow problems, possibly including the contribution of surface tension at the material interfaces. Conceptualization, F.B., L.A.B. and F.C.M. ; Investigation, F.C.M. ; Methodology, F.B., A.C. and F.C.M. ; Supervision, F.B., L.A.B. and A.C.; Validation, F.C.M. ; Writing--original draft, F.C.M. ; Writing--review & editing, F.B., L.A.B., A.C. and F.C.M. All authors have read and agreed to the published version of the manuscript. This research received no external funding. Not applicable. Data availability StatementData will be made available from the corresponding author upon reasonable request. This research is in partnership with the International Research Training Group (IRTG): \"DROPIet Interaction Technologies\" (DROPIT) grant no. GRK 2160/2. The authors declare no conflict of interest. The following abbreviations are used in this manuscript: