Turbulent thermonuclear burning is studied on scales relevant to the explosion of Type Ia supernovae. A scaling law is formulated for turbulent burning in a uniform gravitational field. The steady state turbulent flame speed is Dδt = f(α) √gL in the regime where the Froude number F = D2l/gL ≪ 1; g, L, Dl, and α = ρ0/P1 > 1 are the acceleration, characteristic scale of the problem, normal speed of the laminar flame, and ratio of the densities ahead and behind the flame, respectively; and f ≃ 1 is a universal function. In this regime, the turbulent flame speed does not depend on the laminar speed Dl and on details of burning on scales ≪L. A flame-capturing technique for modeling turbulent burning is described. It is used to numerically study the transition to turbulence and turbulent flame propagation in three dimensions. The results confirm the scaling law. The self-regulating mechanism underlying the scaling law is discussed. In Type Ia supernovae, steady state burning takes place on scales less than the radius of the flame, where the effects of spherical geometry and expansion are small. Larger scales influenced by these effects need to be resolved explicitly. Direct, ab initio three-dimensional numerical simulations of deflagration in supernovae thus become feasible. Effects of spherical geometry and expansion of matter on the propagation of turbulent flames are discussed. The expansion decreases large-scale turbulent motions and reduces the bulk rate of deflagration in a massive carbon-oxygen white dwarf. Results of a large-scale three-dimensional simulation of the deflagration explosion of a Type Ia supernova are presented.