The simulation of rarefied gas flow based on the Boltzmann equation is challenging, especially when the gas mixtures have disparate molecular masses. In this paper, a computationally tractable kinetic model is proposed for monatomic gas mixtures, to mimic the Boltzmann collision operator as closely as possible. The intra- and inter-collisions are modelled separately using relaxation approximations, to correctly recover the relaxation time scales that could span several orders of magnitude. The proposed kinetic model preserves the accuracy of the Boltzmann equation in the continuum regime by recovering four critical transport properties of a gas mixture: the shear viscosity, the thermal conductivity, the coefficients of diffusion and the thermal diffusion. While in the rarefied flow regimes, the kinetic model is found to be accurate when comparing its solutions with those from the direct simulation Monte Carlo method in several representative cases (e.g. one-dimensional normal shock wave, Fourier flow and Couette flow, two-dimensional supersonic flow passing a cylinder and nozzle flow into a vacuum), for binary mixtures with a wide range of mass ratios, species concentrations and different intermolecular potentials. Pronounced separations in species properties have been observed, and the flow characteristics of gas mixtures in shock waves are found to change as the molecular mass ratio increases from 10 to 1000.