  Simulating the flow of two fluid phases in porous media is a challenging
task, especially when fractures are included in the simulation. Fractures may
have highly heterogeneous properties compared to the surrounding rock matrix,
significantly affecting fluid flow, and at the same time hydraulic aperture
that are much smaller than any other characteristic sizes in the domain.
Generally, flow simulators face difficulties with counter-current flow,
generated by gravity and pressure gradients, which hinders the convergence of
non-linear solvers (Newton).
  In this work, we model the fracture geometry with a mixed-dimensional
discrete fracture network, thus lightening the computational burden associated
to an equi-dimensional representation. We address the issue of counter-current
flows with appropriate spatial discretization of the advective fluid fluxes,
with the aim of improving the convergence speed of the non-linear solver. In
particular, the extension of the hybrid upwinding to the mixed-dimensional
framework, with the use of a phase potential upstreaming at the interfaces of
subdomains.
  We test the method across several cases with different flow regimes and
fracture network geometry. Results show robustness of the chosen discretization
and a consistent improvements, in terms of Newton iterations, compared to use
the phase potential upstreaming everywhere.
