In this study, the combined application of the geodesic kernel and Gaussian process regression was explored to estimate nonstationary hydraulic conductivity fields in two-dimensional hydrogeological systems. Specifically, a semianalytical form of the geodesic distance based on the intrinsic geometry of the manifold was derived and used to define positive definite geodesic covariance matrices that were used for the Gaussian process regression. Furthermore, the proposed approach was applied to a series of synthetic hydraulic conductivity estimation problems. The results show that the incorporation of secondary information, such as interpretations of geophysical explorations or geological surveys, can considerably improve the accuracy of the estimation, especially in nonstationary fields. In addition, groundwater flow and solute transport simulations based on estimated hydraulic conductivity fields reveal that the accuracy of the simulation is strongly affected by the inclusion of secondary information. These results suggest that incorporating secondary information into manifold geometry can remarkably improve the accuracy of the estimation and provide new insights into the underlying structure of geological data. This proposed approach has critical implications for hydrogeological applications, such as groundwater resource management, safety assessments, and risk management strategies related to groundwater contamination.