The dynamics of turbulence-driven E x B zonal flows has been systematically studied in fully 3-dimensional gyrokinetic simulations of microturbulence in magnetically confined toroidal plasmas using recently available massively parallel computers. Linear flow damping simulations exhibit an asymptotic residual flow in agreement with recent analytic calculations. Nonlinear global simulations of instabilities driven by temperature gradients in the ion component of the plasma provide key first principles results supporting the physics picture that turbulence-driven fluctuating E x B zonal flows can significantly reduce turbulent transport.