Accurate model predictions of water flow and solute transport in unsaturated soils require a correct representation of relevant mechanisms in a mathematical model, as well as correct solutions of the mathematical equations. Because of the complexity of boundary conditions and the nonlinearity of the processes considered, general solutions of the mathematical equations rely on numerical approximations. We evaluate a number of numerical models (WAVE, HYDRUS_1D, SWAP, MARTHE, and MACRO) that use different numerical methods to solve the flow and transport equations. Our purpose is to give an overview of analytical solutions that can be found for simple initial and boundary conditions and to define benchmark scenarios to check the accuracy of numerical solutions. Included are analytical solutions for coupled transport equations that describe flow and transport in dual-velocity media. The relevance of deviations observed in the analytical benchmarks for more realistic boundary conditions is illustrated using an intercode comparison for natural boundary conditions. For the water flow scenarios, the largest deviations between numerical models and analytical solutions were observed for the case of soil limited evaporation. The intercode differences could be attributed to the implementation of the evaporation boundary condition: the spatial discretization and the internode averaging of the hydraulic conductivity in the surface grid layer. For solute transport, accurate modeling of solute dispersion poses the most problems. Nonlinear and nonequilibrium sorption and coupled transport in pore domains with different advection velocities are in general accurately simulated.