  We consider conservation laws with discontinuous flux where the initial
datum, the flux function, and the discontinuous spatial dependency coefficient
are subject to randomness. We establish a notion of random adapted entropy
solutions to these equations and prove well-posedness provided that the spatial
dependency coefficient is piecewise constant with finitely many
discontinuities. In particular, the setting under consideration allows the flux
to change across finitely many points in space whose positions are uncertain.
We propose a single- and multilevel Monte Carlo method based on a finite volume
approximation for each sample. Our analysis includes convergence rate estimates
of the resulting Monte Carlo and multilevel Monte Carlo finite volume methods
as well as error versus work rates showing that the multilevel variant
outperforms the single-level method in terms of efficiency. We present
numerical experiments motivated by two-phase reservoir simulations for
reservoirs with varying geological properties.
