Fluidized bed reactors are widely used in metallurgical, chemical, food processing, pharmaceutical and energy industries, and wastewater treatment plants due to high heat and mass transfer rates. In practice, these fluidized beds have a wide particle size distribution. Even though uniform-sized particles are charged, they may undergo a reduction in size due to mechanical disintegration, thermal disintegration, and chemical reactions, leading to elutriation of fines in columnar beds. So, it is believed that the tapered fluidized beds are best suited for practical applications as they have a vertical velocity gradient due to an increase in cross-section. It allows fine particles to stay up in the bed, whereas coarse particles stay at the bottom of the bed. Further, it shows different hydrodynamic behavior, especially bubble dynamics, compared to the columnar beds. However, there are very limited studies reported in the literature on bubble behavior, which governs the heat transfer, mass transfer, reaction kinetics, and mixing-segregation phenomena in a tapered fluidized bed. Therefore, the present study aims to investigate the bubble behavior in pseudo-2D gas-solid tapered fluidized beds through physical and mathematical modeling approaches. The fluidization phenomenon was captured using a high-speed camera, and further, the images were analyzed to quantify various parameters. The bubble-based Energy Minimization Multi-Scale (EMMS) drag models have been used to capture the heterogeneity in the bed's flow structure caused by bubbles, as these are better than the homogeneous models. However, all the existing EMMS models contain the implicit interdependency between the bubble size and its drag coefficient. Therefore, a modified bubble-based EMMS drag law is proposed and coupled with the two-fluid Kinetic Theory of Granular Flow (KTGF) model to represent gas-solid flow and solved using Computational Fluid Dynamics (CFD) technique. The effects of taper angle, gas velocity, and particle size on local bubble properties such as its size, shape, rise velocity, and bulk bubble properties such as bubble fraction, bed expansion ratio, and unfluidized area fraction are studied. The bubble coalescence and breakup behavior have also been investigated. Empirical correlations are developed for predicting time-averaged bubble fraction, bed expansion ratio, and unfluidized area fraction. An operating regime map is proposed based on the fraction of unfluidized area. The simulation results obtained using modified EMMS law have been compared with Gidaspow law and present experiments. In comparison with experimental findings, it is found that modified EMMS law is more accurate than the Gidaspow law in predicting the bulk properties.