Turbulent dynamical systems characterized by both a high-dimensional phase space and a large number of instabilities are ubiquitous among many complex systems in science and engineering including climate, material, and neural science. The existence of a strange attractor in the turbulent systems containing a large number of positive Lyapunov exponents results in a rapid growth of small uncertainties from imperfect modeling equations or perturbations in initial values, requiring naturally a probabilistic characterization for the evolution of the turbulent system. Uncertainty quantification (UQ) in turbulent dynamical systems is a grand challenge where the goal is to obtain statistical estimates such as the change in mean and variance for key physical quantities in their nonlinear responses to changes in external forcing parameters or uncertain initial data. In the development of a proper UQ scheme for systems of high or infinite dimensionality with instabilities, significant model errors compared with the true natural signal are always unavoidable due to both the imperfect understanding of the underlying physical processes and the limited computational resources available through direct Monte-Carlo integration. One central issue in contemporary research is the development of a systematic methodology that can recover the crucial features of the natural system in statistical equilibrium (model fidelity) and improve the imperfect model prediction skill in response to various external perturbations (model sensitivity).
   Here we discuss a general mathematical framework to construct statistically accurate reduced-order models that have skill in capturing the statistical variability in the principal directions with largest energy of a general class of damped and forced complex turbulent dynamical systems. There are three stages in the modeling strategy, imperfect model selection; calibration of the imperfect model in a training phase using only data in the more complex perfect model statistics; and prediction of the responses with UQ to a wide class of forcing and perturbation scenarios. The methods are developed under a universal class of turbulent dynamical systems with quadratic nonlinearity that is representative in many applications in applied mathematics and engineering. Several mathematical ideas will be introduced to improve the prediction skill of the imperfect reduced-order models. Most importantly, empirical information theory and statistical linear response theory are applied in the training phase for calibrating model errors to achieve optimal imperfect model parameters; and total statistical energy dynamics are introduced to improve the model sensitivity in the prediction phase especially when strong external perturbations are exerted. The validity of general framework of reduced-order models is demonstrated on instructive stochastic triad models. Recent applications to two-layer baroclinic turbulence in the atmosphere and ocean with combinations of turbulent jets and vortices are also surveyed. The uncertainty quantification and statistical response for these complex models are accurately captured by the reduced-order models with only 2×10² modes in a highly turbulent system with 1×10⁵ degrees of freedom. Less than 0.15% of the total spectral modes are needed in the reduced-order models.