We examine the accuracy of spatial derivatives computed from numerical simulations of supersonic turbulence. Two sets of simulations, carried out using a finite-volume code that evolves the hydrodynamic equations with an approximate Riemann solver and a finite-difference code that solves the Navier-Stokes equations, are tested against a number of criteria based on the continuity equation, including exact results at statistically steady state. We find that the spatial derivatives in the Navier-Stokes runs are accurate and satisfy all the criteria. In particular, they satisfy our exact results that the conditional mean velocity divergence, $langle abla cdot {bs u}|srangle$, where $s$ is the logarithm of density, and the conditional mean of the advection of $s$, $langle {bs u} cdot abla s|srangle$, vanish at steady state for all density values, $s$. On the other hand, the Riemann solver simulations fail all the tests that require accurate evaluation of spatial derivatives, resulting in apparent violation of the continuity equation, even if the solver enforces mass conservation. In particular, analysis of the Riemann simulations may lead to the incorrect conclusion that the $pdv$ work tends to preferentially convert kinetic energy into thermal energy, inconsistent with the exact result that the energy exchange by $pdv$ work is symmetric in barotropic supersonic turbulence at steady state. The inaccuracy of spatial derivatives is a general problem in the post-processing of simulations of supersonic turbulence with Riemann solvers. Solutions from such simulations must be used with caution in post-processing studies concerning the spatial gradients.