Nodal integral methods (NIM) are a class of efficient coarse mesh methods that use transverse averaging to reduce the governing partial differential equation(s) (PDE) into a set of ordinary differential equations (ODE). These ODEs or their approximations are analytically solved, and the solutions are used to develop the set of discrete equations. Since this method depends on transverse averaging, the standard application of NIM gets restricted to domains with boundaries parallel to one of the coordinate axes (2D) or coordinate planes (3D). Hybrid nodal-integral/finite-element method (NI-FEM) has recently been developed to extend the application of NIM to arbitrary domains. NI-FEM is based on the idea that the bulk of the domain and the regions with boundaries parallel to the coordinate axes (2D) or coordinate planes (3D) are discretized using coarse NIM cells (NIM subdomains), and the rest of the domain is discretized using FEM elements (FEM subdomains). The crux of the hybrid NI-FEM is in developing interfacial conditions at the common interfaces between the NIM and the FEM subdomains. Since the discrete variables in the two numerical approaches are different, this requires special treatment of the discrete quantities on the interfaces. We here report the development of hybrid NI-FEM for the 2D, incompressible Navier-Stokes equations (NSE) that are coupled to the energy equation via the Boussinesq approximation. The hybrid scheme is implemented in a parallel framework in Fortran using PETSc. The non-linear system of equations is solved using the Jacobian-free Newton-Krylov (JFNK) approach. The approach is compared with standard FEM to study its efficiency.