AbstractThe use of orthonormal polynomial bases has been found to be efficient in preventing ill‐conditioning of the system matrix in the primal formulation of Virtual Element Methods (VEM) for high values of polynomial degree and in presence of badly‐shaped polygons. However, we show that using the natural extension of a orthogonal polynomial basis built for the primal formulation is not sufficient to cure ill‐conditioning in the mixed case. Thus, in the present work, we introduce an orthogonal vector‐polynomial basis which is built ad hoc for being used in the mixed formulation of VEM and which leads to very high‐quality solution in each tested case. Furthermore, a numerical experiment related to simulations in Discrete Fracture Networks (DFN), which are often characterised by very badly‐shaped elements, is proposed to validate our procedures.