The deterministic approach to modeling groundwater flow overlooks the inherent heterogeneity and uncertainties in aquifer properties and hydrologic conditions. Stochastic models, accounting for spatial variability and uncertainty, offer a more comprehensive solution, vital for informed decision-making in engineering applications, especially at watershed scales. Progress in understanding stochastic unconfined flow has been slower due to the complexity of the problem and the limitations of existing methods. The commonly used stochastic methods such as Monte Carlo and perturbation methods have drawbacks including computational intensity and low parameter variance, especially in complex systems. Addressing these limitations, Kavvas (2003) introduced the 'Master Key' equation, which upscales stochastic dynamical processes to field scale equations. The master key equation converts the governing equations of a stochastic process to their Lagrangian-Eulerian extension of the Fokker-Planck equation (LEFPE) to exact second order. The resulting LEFPE is a linear deterministic partial differential equation that describes the spatio-temporal evolution of the probability density function of state variables. The advantages of this method over traditional stochastic methods are that it provides the ensemble behavior in a single simulation, and it can accommodate large parameter variances. Therefore, this study aims to upscale the two-dimensional stochastic flow equation in a heterogeneous, anisotropic unconfined aquifer with asymmetric properties, accounting for uncertainties in hydraulic conductivity, source/ sink term, and boundary conditions, using the aforementioned LEFPE approach. One-dimensional and two-dimensional numerical models of unconfined aquifer flow were developed using the LEFPE method. The results were compared against Monte Carlo simulations to assess the ensemble average behavior, standard deviation, and probability density functions, thereby providing a comprehensive analysis of the probabilistic description of the stochastic unconfined aquifer flow process modeled by the methodology developed in this study. These results highlight the effectiveness of the LEFPE method in efficiently capturing the evolutionary probabilistic behavior of hydraulic head in a heterogeneous, asymmetric unconfined aquifer.