Review of numerical solution of Richardson-Richards equation for variably saturated flow in soils

Modeling variably saturated flow in the vadose zone is of vital importance to many scientific fields and engineering applications. Richardson–Richards equation (RRE, which is conventionally known as Richards' equation) is often chosen to physically represent the fluxes in the vadose zone when the accurate characterization of the soil water dynamics is required. Being a highly nonlinear partial differential equation, RRE is often solved numerically. Although there are mature software and codes available for simulating variably saturated flow by solving RRE, the numerical solution of RRE is nevertheless computationally expensive. Moreover, sometimes the robustness and the efficiency of RRE‐based models can deteriorate rapidly when certain unfavorable conditions are met. These demerits of RRE hinder its application on large‐scale vadose zone hydrology problems and uncertainty quantification, both of which requires many runs of the prediction model. To address these challenges, the accuracy, convergence, and efficiency of the numerical schemes of RRE should be further improved by testing a wide variety of cases covering different initial conditions, boundary conditions, and soil types. We reviewed and highlighted several critical issues related to the numerical modeling of RRE, including spatial and temporal discretization, the different forms of RREs, iterative and noniterative schemes, benchmark solutions, and available software and codes. Based on the review, we summarize the challenges and future work for solving RRE numerically.