This paper develops a Bayesian optimal experimental design for robot kinematic calibration on ${\mathbb{S}^3 \!\times\! \mathbb{R}^3}$. Our method builds upon a Gaussian process approach that incorporates a geometry-aware kernel based on Riemannian Matérn kernels over ${\mathbb{S}^3}$. To learn the forward kinematics errors via Bayesian optimization with a Gaussian process, we define a geodesic distance-based objective function. Pointwise values of this function are sampled via noisy measurements taken using fiducial markers on the end-effector using a camera and computed pose with the nominal kinematics. The corrected Denavit-Hartenberg parameters are obtained using an efficient quadratic program that operates on the collected data sets. The effectiveness of the proposed method is demonstrated via simulations and calibration experiments on NASA's ocean world lander autonomy testbed (OWLAT).