We present the first detailed three-dimensional ( 3D ) hydrodynamic implicit large eddy simulations of turbulent convection of carbon burning in massive stars . Simulations begin with radial profiles mapped from a carbon burning shell within a 15 \textrm { M } _ { \odot } one-dimensional stellar evolution model . We consider models with 128 ^ { 3 } , 256 ^ { 3 } , 512 ^ { 3 } and 1024 ^ { 3 } zones . The turbulent flow properties of these carbon burning simulations are very similar to the oxygen burning case . We performed a mean field analysis of the kinetic energy budgets within the Reynolds-averaged Navier-Stokes framework . For the upper convective boundary region , we find that the numerical dissipation is insensitive to resolution for linear mesh resolutions above 512 grid points . For the stiffer , more stratified lower boundary , our highest resolution model still shows signs of decreasing sub-grid dissipation suggesting it is not yet numerically converged . We find that the widths of the upper and lower boundaries are roughly 30 % and 10 % of the local pressure scale heights , respectively . The shape of the boundaries is significantly different from those used in stellar evolution models . As in past oxygen-shell burning simulations , we observe entrainment at both boundaries in our carbon-shell burning simulations . In the large Péclet number regime found in the advanced phases , the entrainment rate is roughly inversely proportional to the bulk Richardson number , Ri _ { B } ( \propto Ri { { } _ { B } } ^ { - \alpha } , 0.5 \lesssim \alpha \lesssim 1.0 ) . We thus suggest the use of Ri _ { B } as a means to take into account the results of 3D hydrodynamics simulations in new 1D prescriptions of convective boundary mixing .