This thesis presents high-order, discontinuous Galerkin (DG) discretizations of the Reynolds-Averaged Navier-Stokes (RANS) equations and an output-based error estimation and mesh adaptation algorithm for these discretizations. In particular, DG discretizations of the RANS equations with the Spalart-Allmaras (SA) turbulence model are examined. The dual consistency of multiple DG discretizations of the RANS-SA system is analyzed. The approach of simply weighting gradient dependent source terms by a test function and integrating is shown to be dual inconsistent. A dual consistency correction for this discretization is derived. The analysis also demonstrates that discretizations based on the popular mixed formulation, where dependence on the state gradient is handled by introducing additional state variables, are generally asymptotically dual consistent. Numerical results are presented to confirm the results of the analysis. The output error estimation and output-based adaptation algorithms used here are extensions of methods previously developed in the finite volume and finite element communities. In particular, the methods are extended for application on the curved, highly anisotropic meshes required for boundary conforming, high-order RANS simulations. Two methods for generating such curved meshes are demonstrated. One relies on a user-defined global mapping of the physical domain to a straight meshing domain. The other uses a linear elasticity node movement scheme to add curvature to an initially linear mesh. Finally, to improve the robustness of the adaptation process, an "unsteady" algorithm, where the mesh is adapted at each time step, is presented. The goal of the unsteady procedure is to allow mesh adaptation prior to converging a steady state solution, not to obtain a time-accurate solution of an unsteady problem. Thus, an estimate of the error due to spatial discretization in the output of interest averaged over the current time step is developed. This error estimate is then used to drive an h-adaptation algorithm. Adaptation results demonstrate that the high-order discretizations are more efficient than the second-order method in terms of degrees of freedom required to achieve a desired error tolerance. Furthermore, using the unsteady adaptation process, adaptive RANS simulations may be started from extremely coarse meshes, significantly decreasing the mesh generation burden to the user.