This manuscript presents a method to identify the occurrence of Equatorial Plasma Bubbles (EPBs) with data gathered from receivers of Global Navigation Satellite System (GNSS). This method adapts a previously existing technique to detect Medium Scale Travelling Ionospheric Disturbances (MSTIDs), which focus on the 2nd time derivatives of total electron content estimated from GNSS signals (2DTEC). Results from this tool made possible to develop a comprehensive analysis of the characteristics of EPBs. Analyses of the probability of occurrence, effective time duration, depth of the depletion and total disturbance of the EPBs show their dependence on local time and season of the year at global scale within the latitude belt from 35°N to 35°S for the descending phase of solar cycle 23 and ascending phase of solar cycle 24, 2002–2014. These results made possible to build an EPBs model, bounded with the Solar Flux index, that simulates the probability of the number of EPBs and their characteristics expected for a representative day at given season and local time (LT). The model results provided insight into different important aspects: the maximum occurrence of bubbles take place near the equatorial anomaly crests, asymmetry between hemispheres and preferred longitudes with enhanced EPBs activity. Model output comparisons with independent observations confirmed its soundness.