Stellarators are inherently steady-state and disruption-free. However, three-dimensional geometry significantly increases the complexity of both equilibrium calculation and configuration optimization. While global three-dimensional magnetohydrodynamic equilibria are generally solved numerically, asymptotic expansions provide a powerful tool for studying stellarator configurations. Despite successes with near-axis expansion, the expansion near an arbitrary flux surface presents unresolved theoretical challenges. We investigate the vacuum field constraints in Boozer coordinates and develop a near-surface expansion framework. This approach is validated on a rotating elliptical boundary stellarator and a precise quasiaxisymmetric stellarator. The results from the first-order expansion agree with the global solutions in high precision, demonstrating a quadratic convergence rate as theoretically predicted. Using a recursive first-order expansion, we can derive approximate global solutions that are considerably far away from the reference surface but remain consistent with the 3D MHD equilibrium code. The shape of the magnetic axis can also be quickly estimated with a few recursive calculations. This work generalizes the asymptotic expansion method to an arbitrary flux surface and provides a new perspective for stellarator configuration studies.