  The outer envelopes of massive ($M\gtrsim10\,M_{\odot}$) stars exhibit large
increases in opacities from forests of lines and ionization transitions
(particularly from iron and helium) that trigger near-surface convection zones.
One-dimensional models predict density inversions and supersonic motions that
must be resolved with computationally intensive 3D radiation hydrodynamic (RHD)
modeling. Only in the last decade have computational tools advanced to the
point where ab initio 3D models of these turbulent envelopes can be calculated,
enabling us to present five 3D RHD Athena++ models (four previously published
and one new 13$M_{\odot}$ model). When convective motions are sub-sonic, we
find excellent agreement between 3D and 1D velocity magnitudes, stellar
structure, and photospheric quantities. However when convective velocities
approach the sound speed, hydrostatic balance fails as the turbulent pressure
can account for 80% of the force balance. As predicted by Henyey, we show that
this additional pressure support leads to a modified temperature gradient which
reduces the superadiabaticity where convection is occurring. In addition, all
five models display significant overshooting from the convection in the Fe
convection zone. As a result, the turbulent velocities at the surface are
indicative of those in the Fe zone. There are no confined convection zones as
seen in 1D models. In particular, helium convection zones seen in 1D models are
significantly modified. Stochastic low frequency brightness variability is also
present in the 13$M_{\odot}$ model with comparable amplitude and characteristic
frequency to observed stars.
