The Wavelet Adaptive Multiresolution Representation (WAMR) method provides a robust method for controlling spatial grid adaption --- fine grid spacing in regions of a solution requiring high resolution (i.e. near steep gradients, singularities, or near-singularities) and using much coarser grid spacing where the solution is slowly varying. The sparse grids produced using the WAMR method exhibit very high compression ratios compared to uniform grids of equivalent resolution. Subsequently, a wide range of spatial scales often occuring in continuum physics models can be efficiently captured. Furthermore, the wavelet transform provides a direct measure of local error at each grid point, effectively producing automatically verified solutions. In this work, the WAMR algorithm is parallelized using an MPI-based domain decomposition approach suitable for a wide range of distributed-memory parallel architectures. The method is applied to the solution of the unsteady, compressible, reactive Navier-Stokes equations and includes detailed diffusive transport and chemical kinetics models. Performance of the method is examined on several test problems in 1-, 2-, and 3-dimensions. The WAMR method shows an impressive compression of the solution, reducing the number of collocation points used by factors of up to $O(10,000)$. Excellent parallel performance is also seen, with near-linear scaling for up to 512 processor cores tested. The WAMR method is then applied to three high-speed compressible flow problems. The first is the shock and reshock of a thin, dense, gas layer in air, and the subsequent Richtmyer-Meshkov instability which develops. The second problem is a cellular detonation in a hydrogen-oxygen-argon mixture. Lastly, the ignition and combustion of a hydrogen bubble by a shockwave in air is simulated. In all cases, results agree favourably with experimental and previous computational results.