We report on a few recent results related to thermal convection in a fluid layer overlying a saturated porous media based on the Navier-Stokes-Darcy-Boussinesq (NSDB) model with appropriate interface boundary conditions. The existence of global in time weak solution for the NSDB system together with a weak-strong uniqueness result are presented first. The stability of the pure conduction state at small Rayleigh number is introduced next. The loss of stability of the pure conduction state as the Rayleigh number crosses a threshold value is studied via a hybrid approach that combines analysis with numerical computation. In particular, we discover that the transition between shallow and deep convection associated with the variation of the ratio of free-flow to porous media depth is accompanied by the change of the most unstable mode from the lowest possible horizontal wave number to higher wave numbers, which could occur with variation of the height ratio as well as the Darcy number and the ratio of thermal diffusivity among others. Numerical methods that decouples the heat, the free flow, the porous media flow while maintaining energy stability are presented as well.