Data-driven reduced-order models often fail to accurately forecast high-dimensional nonlinear systems that are sensitive along coordinates with low-variance. Such systems are encountered frequently in shear-dominated fluid flows where non-normality plays a significant role in the growth of disturbances. In order to address these issues, we employ ideas from active subspaces to find low-dimensional coordinates for model reduction that balance adjoint-based information about the system’s sensitivity with the variance of states along trajectories. The resulting method, which we refer to as covariance balancing reduction using adjoint snapshots (CoBRAS), is identical to balanced truncation with state and adjoint-based gradient covariance matrices replacing the system Gramians and obeying the same transformation laws. Here, the extracted coordinates are associated with an oblique projection that can be used to construct Petrov-Galerkin reduced-order models. We also observe that the reduced coordinates can be computed relying on inner products of state and gradient samples alone, allowing us to find rich nonlinear coordinates by replacing the inner product with a kernel function. We demonstrate these techniques on an axisymmetric jet flow simulation with 100,000 state variables.