Abstract: Increasing the polynomial order of basis functions is an effective strategy for enhancing the accuracy of finite element-based numerical simulations in geophysics. A conventional finite element method typically employs the basis functions with a uniform polynomial order across the entire computational domain. As a result, a higher order significantly increases computational cost. To address this challenge, a finite element method based on hierarchical basis functions is proposed. This method enables local p-refinement specifically for surface elements by using hierarchical basis functions and thereby improves the accuracy of numerical forward modeling without compromising computational efficiency. The validation through audiomagnetotelluric forward modeling on three typical geoelectric models demonstrates that the proposed method significantly enhances the accuracy of forward modeling with only a small increase in computation time, offering a novel approach for efficient numerical electromagnetic forward modeling.