• For the first time, we simulate the bursty bulk flow events with kinetic physics embedded in a global magnetosphere model. • The electron velocity distributions demonstrate different anisotropy features at different locations surrounding the BBF. • Energy dependent electron pitch angle distribution evolutions are identified in the simulation.