Highly variable flux surface averaged heat fluxes are resolved in gyrokinetic simulations of ion temperature gradient (ITG) turbulence, even in large systems. Radially propagating fronts or avalanches are also seen. Their propagation lengths in gyroradii and relative amplitude remain constant as simulation size is increased, so the avalanches appear to result from local dynamics, rather than global relaxation events. For the Cyclone [Dimits , Phys. Plasmas 7, 969 (2000)] case, the avalanche propagation direction is found to depend on the sign of the shearing rate. A mechanism for avalanche propagation based on the advection of turbulence tilted by the shear flows is proposed: The Cyclone linear ITG dispersion relation explains the propagation direction of tilted vortices. It also explains why there is no such preferred direction in a simulation with reduced magnetic shear. The paper explores several models for these bursts. First, certain types of models based on nonlinear heat diffusion equations are ruled out. A different type of one-dimensional (1D) model, introduced in Benkadda [Nucl. Fusion 41, 995 (2001)], yields much better qualitative and quantitative agreement. However, the 1D model cannot explain the directionality of the bursts, even though it includes the features typically considered important for burst propagation. A symmetry-breaking term is necessary. An additional term is included to reproduce the wave dispersion with respect to radial wavenumber, and this is shown to be sufficient to reproduce the favored direction for burst propagation.