Accurate modeling of variably saturated flow (VSF) in fractured porous media with the discrete fracture-matrix (DFM) model is a computationally challenging problem. The applicability of DFM model to VSF in real field studies at large space and time scales is often limited, not only because it requires detailed fracture characterization, but also as it involves excessive computational efforts. We develop an efficient numerical scheme to solve the Richards equation in discretely fractured porous media. This scheme combines the mixed hybrid finite element method for space discretization with the method of lines for time integration. The fractures are modeled as lower-dimensional interfaces (1D), within the 2D porous matrix. We develop a new mass-lumping (ML) technique for the fractures to eliminate unphysical oscillations and convergence issues in the solution, which significantly improves efficiency, enabling larger field applications. The proposed new scheme is validated against a commercial simulator for problems involving water table recharge at the laboratory scale. The computational efficiency of the developed scheme is examined on a challenging problem for water infiltration in fractured dry soil, and compared with standard numerical techniques. We show that the ML technique is crucial to improve robustness and efficiency, which outperforms the commonly used methods that we tested. The applicability of our method is then demonstrated in a study concerning the effect of climate change on groundwater resources in a karst aquifer/spring system in El Assal, Lebanon. Simulations, including recharge predictions under climate change scenarios, are carried out for about 80 years, up to 2099. This study demonstrates the applicability of our proposed scheme to deal with real field cases involving large time and space scales with high variable recharge. Our results indicate that the water-table level is sensitive to the presence of fractures, where neglecting fractures leads to an overestimation of the available groundwater amount. The proposed numerical approach is generic for DFM model and can be extended to different 2D and 3D finite-element frameworks.