We describe a novel and general framework for solving advection-diffusion equations using finite volume weighted essentially non oscillatory (WENO) techniques on general computational meshes. Such techniques are able to handle advective and (degenerate) diffusive behavior, even when the solution develops shocks or steep fronts. We discuss a robust procedure for producing accurate stencil polynomial approximations and a recently developed multilevel WENO (ML-WENO) reconstruction. It combines stencil polynomials of various degrees (e.g., more than two degrees and including constant polynomials) defined on any set of stencils (e.g., not hierarchically arranged). The nonlinear weighting biases the reconstruction away from both inaccurate oscillatory polynomials crossing a shock or steep front and smooth polynomials of low degree, thereby selecting the smooth polynomial(s) of maximal degree of approximation. We apply these ideas to develop a preliminary finite volume scheme for solving the Richards equation, which models unsaturated flow in porous media. Numerical tests of rainwater infiltration show the advantages of using higher order finite volume methods and multilevel WENO reconstructions.