The coupling of free incompressible flow with porous media flow has recently garnered tremendous interest due to its wide range of applications in science, engineering, and industry. It can be used to model hydro-geological mechanics, soil pollution simulation, bio-hydrodynamics, oil drilling and production engineering, industrial filtration, etc. The fluid flow is modeled through Navier-Stokes equations and Darcy’s law governs the porous media flow, and is coupled through certain transmission conditions at the interface. The project aims to investigate this coupled Navier-Stokes-Darcy flow problem with the Beaver-Joseph-Staffman interface condition, which presents numerous mathematical and numerical challenges. As a result, the implementation of the coupled nonlinear discrete problem for a small mesh size is complicated, and the computational cost is significantly higher. Our objective is to analyze numerical schemes, resulting in linearized decoupled models which can be solved separately, leading to easy and efficient implementations. In the proposed problem, two distinct flows are present in the interface region, which causes the system to use either a diffuse or sharp interface to model these free boundaries. However, numerically handling changes in the interface poses a difficult task, as the evolution process may involve the disintegration or merging of interfaces. As the solution varies rapidly over this interface region with time, capturing these free boundary changes is crucial; otherwise, one may observe spurious numerical approximations. Therefore, to capture these changes at the interface, we plan to apply DG methods. These methods are stable, high-order accurate, and locally conservative, even in convection-dominated regimes, and easily handle meshes with hanging nodes, elements of various types and shapes, and local spaces of different orders. Additionally, we propose a combination of two-grid and DG methods for a more cost-efficient approach to solving the proposed nonlinear problem. This approach entails solving the entire nonlinear problem on a coarse mesh, followed by solving a simplified linear version on a finer mesh, ultimately reducing the problem’s size that needs to be solved. The incompressibility condition and mathematical equations pose significant challenges in finite element analysis. To address this, we focus on applying DG stabilized penalty methods to our model, which have proven effective in previous studies. To validate these schemes, we will analyze their formulations for well-posedness and stability, conduct error analysis, and perform numerical computations. By doing so, we hope to understand better the proposed model and the various finite element schemes available. The outcomes of our study will aid in filling critical gaps in the theoretical verification of numerical and experimental results pertaining to the model and other similar kinds of coupled Navier-Stokes models that lack relevant estimates.