A new multifidelity method is developed for nonlinear orbit uncertainty propagation. This approach guarantees improved computational efficiency and limited accuracy losses compared with fully high-fidelity counterparts. The initial uncertainty is modeled as a weighted sum of Gaussian distributions whose number is adapted online to satisfy the required accuracy. As needed, univariate splitting libraries are used to split the mixture components along the direction of maximum nonlinearity. Differential Algebraic techniques are used to propagate these Gaussian kernels and compute a measure of nonlinearity required for the split decision and direction identification. Taylor expansions of the flow of the dynamics are computed using a low-fidelity...