As a simplified model for subsurface flows elliptic equations may be utilized. Insufficient measurements or uncertainty in those are commonly modeled by a random coefficient, which then accounts for the uncertain permeability of a given medium. As an extension of this methodology to flows in heterogeneous\fractured\porous media, we incorporate jumps in the diffusion coefficient. These discontinuities then represent transitions in the media. More precisely, we consider a second order elliptic problem where the random coefficient is given by the sum of a (continuous) Gaussian random field and a (discontinuous) jump part. To estimate moments of the solution to the resulting random partial differential equation, we use a pathwise numerical approximation combined with multilevel Monte Carlo sampling. In order to account for the discontinuities and improve the convergence of the pathwise approximation, the spatial domain is decomposed with respect to the jump positions in each sample, leading to pathdependent grids. Hence, it is not possible to create a nested sequence of grids which is suitable for each sample path a-priori. We address this issue by an adaptive multilevel algorithm, where the discretization on each level is sample-dependent and fulfills given refinement conditions.