Nonlinear uncertainty propagation is of critical importance in many application fields of astrodynamics. In this paper, a framework combining the differential algebra technique and the Gaussian mixture model method is presented to accurately propagate the state uncertainty of a nonlinear system. A high-order Taylor expansion of the final state with respect to the initial deviations is firstly computed with the differential algebra technique. Then the initial uncertainty is split to a Gaussian mixture model. With the high-order state transition polynomial, each Gaussian mixture element is propagated to the final time, forming the final Gaussian mixture model. Through this framework, the final Gaussian mixture model can include the effects of high-order terms during propagation and capture the non-Gaussianity of the uncertainty, which enables a precise propagation of probability density. Moreover, the manual derivation and integration of the high-order variational equations is avoided, which makes the method versatile. The method can handle both the application of nonlinear analytical maps on any domain of interest and the propagation of initial uncertainties through the numerical integration of ordinary differential equation. The performance of the resulting tool is assessed on some typical orbital dynamic models, including the analytical Keplerian motion, the numerical J 2 perturbed motion, and a nonlinear relative motion.