In this research work, the effects of perturbations due to resonant geopotential harmonics on the commensurate semi-major axis \(a_{c}^{{}}\)of the GPS orbits are being studied and analyzed. A second order approximation of the commensurate semi-major axis due to the higher zonal perturbations (up to \({J_6}\)) is obtained iteratively, starting with the unperturbed case as the zeroth-order. In the resonance cases a certain number of tesseral harmonics can produce effects of large amplitude and very long period, these important tesseral harmonics can be computed by Kaula’s resonant perturbation theory [1]. Using Lagrange’s planetary equations, the rate of change of the semi-major axis is computed. The drift rate as a function of longitudinal position is derived and plotted, revealing the existence of four roots which can be classified as stable positions at \({25.94^ \circ }\)E, \({205.9^ \circ }\)E and metastable positions at \({118.7^ \circ }\)E, \({298.7^ \circ }\)E. Motion about these stable and metastable points is investigated using the Poincare section method. The analysis reflects the existence of periodic, quasi-periodic, and even chaotic orbits in the vicinity of these points.