The 2D Boussinesq system is potentially relevant to the study of atmospheric and oceanographic turbulence, as well as other astrophysical situations where rotation and stratification play a dominant role. In fluid mechanics, the 2D Boussinesq system is commonly used in the field of buoyancy-driven flow. It describes the motion of incompressible inhomogeneous viscous fluid subject to convective heat transfer under the influence of gravitational force. It is well-known that the 2D Boussinesq equations are closely related to 3D Euler or Navier-Stokes equations for incompressible flow, and it shares a similar vortex stretching effect as that in the 3D incompressible flow. In fact, in vortex formulation, the 2D inviscid Boussinesq equations are formally identical to the 3D incompressible Euler equations for axisymmetric swirling flow. Therefore, the qualitative behaviors of the solutions to the two systems are expected to be identical. Better understanding of the 2D Boussinesq system will undoubtedly shed light on the understanding of 3D flows. In this talk, I will discuss some recent results concerning global existence, uniqueness and asymptotic behavior of classical solutions to initial boundary value problems for 2D Boussinesq equations with partial viscosity terms on bounded domains for large initial data.