Stochastic field distortions caused by atmospheric turbulence are a fundamental limitation to the astrometric accuracy of ground-based imaging. This distortion field is measurable at the locations of stars with accurate positions provided by the Gaia DR2 catalog; we develop the use of Gaussian process regression (GPR) to interpolate the distortion field to arbitrary locations in each exposure. We introduce an extension to standard GPR techniques that exploits the knowledge that the 2-dimensional distortion field is curl-free. Applied to several hundred 90-second exposures from the Dark Energy Survey as a testbed, we find that the GPR correction reduces the variance of the turbulent distortions $\approx12\times$, on average, with better performance in denser regions of the Gaia catalog. The RMS per-coordinate distortion in the $riz$ bands is typically $\approx7$ mas before any correction, and $\approx2$ mas after application of the GPR model. The GPR astrometric corrections are validated by the observation that their use reduces, from 10 to 5 mas RMS, the residuals to an orbit fit to $riz$-band observations over 5 years of the $r=18.5$ trans-Neptunian object Eris. We also propose a GPR method, not yet implemented, for simultaneously estimating the turbulence fields and the 5-dimensional stellar solutions in a stack of overlapping exposures, which should yield further turbulence reductions in future deep surveys.