We present an in-depth and systematic analysis of a sample of 20 powerful blazars, including 12 BL Lacs and 8 flat spectrum radio quasars, utilizing Fermi/LAT observations from the period 2008-2018 using various analysis tools such as flux distribution, symmetry analysis, and time series analysis. Results show that blazars with steeper γ-ray spectral indexes are found to be more variable; and the γ-ray flux distribution closely resemble both normal and lognormal probability distribution functions. The statistical variability properties of the sources as studied by power spectral density analysis are consistent with flicker noise (P (ν) ∝ 1/ν) -an indication of long-memory processes at work. Statistical analysis of the distribution of flux rise and decay rates in the light curves of the sources, aimed at distinguishing between particle acceleration and energy dissipation timescales, counter-intuitively suggests that both kinds of rates follow a similar distribution and the derived mean variability timescales are in the order of a few weeks. The corresponding emission region size is used to constrain location of γ-ray production sites in the sources to be a few parsecs. Additionally, using Lomb-Scargle periodogram aided with extensive Monte Carlo simulations, we detected year timescale quasi-periodic oscillations in the sources S5 0716+714, Mrk 421, ON +325, PKS 1424-418 and PKS 2155-304; and the detection significance was computed taking proper account of the red-noise and other artifacts inherent in the observations. We explain our results in the light of current blazar models with relativistic shocks propagating down the jet viewed close to the line of sight.