The dynamics of soil moisture in the vadose zone are mathematically modeled by a nonlinear partial differential equation (PDE), which is referred to as the Richardson-Richards equation (RRE;Richards, 1931;Richardson, 1922). Conventionally, theoretical methods for deriving the RRE are rooted in physics laws, including the conservation of mass and the Buckingham-Darcy law (Buckingham, 1907) with conditional simplification assumptions (Narasimhan, 2007;Pinder & Celia, 2006). A thorough understanding of its physical processes is required to determine the soil moisture flow governing equation. However, due to the complexity and uncertainties in the vadose zone, such as the preferential flow through soil macropores (Beven & Germann, 2013;Nimmo, 2021), soil heterogeneity (Elkateb et al., 2003, and changing boundary conditions (Van Dam & Feddes, 2000), a physical principle-based rigorous derivation is complicated and even impossible. In addition, the applications of the RRE rely on accurate measurements of constitutive relationship model parameters, including the water retention curves and hydraulic conductivity functions (Madi et al., 2018). Nevertheless, the determination of soil hydraulic properties is time-consuming, laborious, and costly (Bitterlich et al., 2004). In addition, the experimental scale mismatch problem can worsen the results (Hopmans et al., 2002).With the continuous development of statistical machine learning algorithms, data-driven methods have achieved great success (Remesan & Mathew, 2015). Recently, a seminal breakthrough involved the identification of physical equations through data-driven approaches. Data-driven approaches aim to directly identify physical equations in ordinary differential equations or PDEs from the observed data. Related research traces back to Bongard and Lipson (2007) and Schmidt and Lipson (2009). They successfully used symbolic regression (Koza & Koza, 1992) to discover the underlying dynamic equations from captured motion-tracking data. Subsequently, more data-driven methods, including sparse regression (Brunton et al., 2016;Tibshirani, 1996) and inverse methods, such as Gaussian process regression (Raissi et al., 2017), data assimilation (Chang & Zhang, 2019a), and physics-informed neural networks (Raissi et al., 2019), have also been successfully applied to physical equation discovery.Considering these novel data-driven approaches have shown promising potential in discovering physics equations from data, we raise the following scientific question. Leveraging the powerful learning capability of data-driven methods, can the actual soil moisture flow governing equation be recovered from data, and can the automatic Abstract Conventionally, soil moisture dynamics are mathematically modeled by the Richardson-Richards equation, whose derivation is based on the conservation of mass and the Buckingham-Darcy law. However, it is complicated and even impossible to finish such rigorous derivations based on physical principles due to the complexity and uncertainties in the vadose zon...