This paper shows an application of a vertex-based Compatible Discrete Operator (CDO) scheme to simulate the Richards equation on meshes using polyhedral cells. Second order spatial accuracy is reached for a verification test case using two meshes made of different types of polyhedral cells. A second validation test case dealing with a variably saturated soil inside a vertical column has been simulated, with three advected passive pollutants being released. Results are in good agreement with the reference for both types of mesh. MPI scalability has also been assessed at scale up to 98,304 cores of a Cray X30 for a 1.7 B cell mesh for the simulation of the aforementioned case, and nearly linear speed-up is observed. Finally the implementation of hybrid MPI-OpenMP has been tested on 2 types of Intel Xeon Phi Knights Landing co-processors, showing that the best configuration for the code is obtained when all the physical cores are filled by MPI tasks and the 4 hyperthreads by OpenMP threads.