The vadose zone is the portion of the subsurface above the water table and its pore space usually contains air and water. Due to the presence of infiltration, erosion, plant growth, microbiota, contaminant transport, aquifer recharge, and discharge to surface water, it is crucial to predict the transport rate of water and other substances within this zone. However, flow in the vadose zone has many complications as the parameters that control it are extremely sensitive to the saturation of the media, leading to a nonlinear problem. This flow is referred as unsaturated flow and is governed by Richards equation. Analytical solutions for this equation exists only for simplified cases, so most practical situations require a numerical solution. Nevertheless, the nonlinear nature of Richards equation introduces challenges that causes numerical solutions for this problem to be computationally expensive and, in some cases, unreliable. High order mimetic finite difference operators are discrete analogs of the continuous differential operators and have been extensively used in the fields of fluid and solid mechanics. In this work, we present a numerical approach involving high order mimetic operators along with a Newton root-finding algorithm for the treatment of the nonlinear component. Fully-implicit time discretization scheme is used to deal with the problem’s stiffness.