Understanding nitrate migration through the deep vadose zone is essential for aquifer vulnerability assessments. The effect of variability of physical properties of the deep vadose zone on nitrate transport has been scarcely explored. Recently, deep nitrate storage profiles were determined in the vadose zone of the Loess Plateau of China. Using these observations along with measured soil properties, this study investigates the effect of loess vertical heterogeneity on water movement and nitrate transport through the deep vadose zone. Models of different complexity were established and calibrated. First, a simple piston flow and nitrate mass balance approach was calibrated to the observed nitrate storage. The results indicate that the total nitrate storage is estimated well, while the estimation of the distribution of nitrate is relatively poor. Subsequently, Richards’ equation and the Advection-Dispersion equation were evaluated. Three different conceptualizations of the numerical models were calibrated against deep vadose zone nitrate and water content observations: (1) one-layer model assuming homogenous loess vadose zone; (2) a model that considers a hydraulic conductivity (Ks) decay function and (3) a model where the Miller-Miller scaling factors are prescribed to account for changes of the hydraulic functions with depth. Accounting for the vertical Ks decay in the numerical models improved water flow performances. The study reveals the adequacy of implementing water flow and nitrate transport numerical models together with a simple representation of the vertical loess variability, for simulating nitrate migration in loess deep vadose zone environments.