Mathematical models and simulations of processes related to the interactions between the root system and the surrounding soil can serve as valuable tools for understanding how plant roots influence soil water distribution and for selecting advantageous root traits. These mathematical models typically represent the root system architecture as a 1D branched network embedded in a 3D soil porous media, assuming that the root diameter is smaller than the main spatial dimension of the soil. The water flow in the vadose zone around the root plant system and in the root xylem can be modelled through the non-linear time-dependent possible degenerate Richards equation. In this talk, we present a new mathematical approach to numerically solve this strongly non-linear coupled problem that exploits the mixed Virtual Element Methods to ensure mass conservation. Moreover, since we consider an evolving root system architecture, we rely on a 3D-1D coupling strategy that eliminates the necessity for mesh conformity at the interface. Consequently, this approach mitigates any additional costs associated with the need for remeshing as the network geometry changes.