Generally, the fluxes \(\mathbf {F}_{i + \frac{1}{2}}\)  and \(\mathbf {F}_{i - \frac{1}{2}}\)  are computed by solving Riemann problems at cell boundary \(x_{i+\frac{1}{2}}\)  with input data \(\mathbf {U}_i\)  and \(\mathbf {U}_{i+1}\)  and at cell boundary \(x_{i-\frac{1}{2}}\)  with input data \(\mathbf {U}_{i-1}\)  and \(\mathbf {U}_{i}\) , respectively.
Here, for fluxes at interfaces between cells with equal phases a standard Riemann solver can be used, see e.g. Toro[1]} for an overview. We specify our choice when presenting the numerical results. For the cell boundaries with phase transition, i.e., the interfaces between cells \(I_{k-1}\)  and \(I_k\)  as well as cells \(I_k\)  and \(I_{k+1}\) , we employ special Riemann solvers. For completeness, we provide the description in the appendix . The focus of this contribution is on the choice of the time stepping scheme, i.e., on which instance of time to evaluate the input arguments \(\mathbf {U}_{i \pm 1}\)  and \(\mathbf {U}_{i}\) .
